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ABSTRACT 


The  three-dimensional  singular  stress  field  near  the  terminal 
point  0  of  the  crack  front  edge  at  the  surface  of  an  elastic  body  is 
investigated.  Displacements  are  assumed  to  be  of  the  form  r^p^F(6,0) 
where  spherical  coordinates  r,  9,  0  are  used,  and  where  p  is  the  distance 
from  the  singularity  line  (crack  front  edge  or  notch  edge)  and  p  is  a 
given  constant.  The  variational  principle  governing  the  displacement 
distribution  on  a  unit  sphere  about  point  0  is  derived  from  the  dif¬ 
ferential  equations  of  equilibrium,  and  more  directly,  from  the  po-^ 
tential  energy.  A  finite  element  method  developed  on  the  unit  sphere 
is  used  to  reduce  the  problem  to  the  form  k(X)x  *  0  where  x  is  the 
column  matrix  of  the  nodal  values  of  the  displacements  on  the  unit 
sphere  and  k(X.)  is  a  square  matrix,  all  coefficients  of  which  are 
quadratic  polynomials  in  X-  It  is  proved  that  the  variational  princi¬ 
ple  as  well  as  the  matrix  k  must  be  nonsymme tr ic ,  which  implies  that 
complex  eigenvalues  \  are  possible.  Several  numerical  and  analytical 
solutions  are  compared  and  agree  closely  with  the  present  work.  By 
energy  flux  arguments  it  is  found  that  the  front  edge  of  a  propagating 
crack  must  terminate  at  the  surface  obliquely  at  a  certain  angle, 
whose  dependence  upon  the  inclination  of  a  crack  plane  is  also  solved. 

The  angle  is  the  same  for  Modes  II  and  III,  but  different  for  Mode  I. 

For  Mode  I,  the  surface  point  trails  behind  the  interior  of  the  pro¬ 
pagating  crack,  while  for  Modes  II  and  III  it  moves  ahead.  Consequently, 
a  combination  of  Mode  I  with  Modes  II  and  III  is  impossible  at  the 
surface  terminal  point  of  a  propagating  crack  whose  plane  is  orthogonal. 


When  the  plane  is  inclined,  the  three  intensify  factors  can  com-  j 

| 

bine  only  in  certain  fixed  ratios.  The  crack  edge  angle  is  a  ; 

function  of  the  angle  of  the  crack  plane.  Some  results  are  also 
presented  for  notches  and  for  cracks  that  intersect  a  two-material 
interface. 
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INTRODUCTION 

Crack  propagation  in  thin  sheets  is  undoubtedly  influenced  by 
the  surface  termination  of  the  crack  front  edge,  where  the  planar 
elasticity  solution  for  the  crack-tip  singularity  does  not  apply  and 
the  singular  stress  field  is  of  three-dimensional  nature.  A  similar 
situation  arises  when  the  crack  front  edge  intersects  a  two-material 
interface.  Knowledge  of  the  three-dimensional  singularity  is  needed 
to  determine  the  curved  shape  of  the  crack  front  edge  across  a  thin 
sheet  or  plate,  and  the  energy  release  rate  for  the  advance  of  the 
crack  front  edge  as  a  whole. 

The  three-dimensional  displacement  field  near  the  terminal  point 
of  the  crack  front  edge  at  the  surface  of  an  elastic  body  is  investi¬ 
gated  in  Chapter  I  using  spherical  corrdinates  r,  9,  0.  The  basic 
idea  of  the  present  method  of  solution  consists  of  an  extension  of  the 
Kne in-William's  method  [1,2]  by  assuming  that  all  three  displacements 
are  of  the  form  r^pPF(9,0)  [3],  where  p  is  the  distance  from  the 
singularity  line,  such  as  the  crack  front  edge  or  notch  edge;  p  is  a 
given  constant;  F(9,0)  is  an  arbitrary  function  of  the  coordinates 
9  and  0;  r  is  the  distance  from  the  terminal  point  0;  and  \  is  the 
strength  singularity  exponent.  Similar  techniques  for  solving  two- 
dimensional  axisymmetric  problems  have  also  been  employed  in  Refs. 

[it- ,5  ,6].  A  partly  similar  approach  has  been  used  by  Swedlow  and 
Karabin  [7]* 

The  variational  equation  governing  the  displacement  distribution 
on  a  unit  sphere  about  the  sirg”1arity  point  0  is  derived  from  the 
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differential  equations  of  equilibrium  and  boundary  conditions,  and 
an  alternate  derivation  is  obtained  more  directly  from  the  potential 
energy.  It  is  proved  that  the  variational  principle  must  be  non- 
symmetric  and  therefore  complex  eigenvalues  X  are  possible.  Thus, 
the  variational  principle  is  of  general  applicability  since  it  can 
handle  crack  and  notches  of  any  orientation  and  size  as  well  as 
problems  of  two  dissimilar  materials. 

The  variational  equation  derived  in  Chapter  I  is  suitable  for 
solution  by  numerical  techniques.  In  Chapter  II,  the  finite 
element  method  is  applied  and  the  problem  is  reduced  to  the  form 
Jc(X)x_  *  0,  where  x  is  the  column  matrix  of  the  nodal  values  of  the 
displacements  on  the  unit  sphere  and  Jc(X)  is  a  banded  square  matrix, 
all  coefficients  of  which  are  quadratic  polynomials  in  X.  The  method 
of  search  for  the  eigenvalue  involves  a  conversion  to  a  non-homoge- 
neous  system  of  equations  and  an  iteration  scheme.  This  method  has 
been  used  in  connection  with  other  problems  which  lead  to  equations 
of  this  type  [5].  Convergence  patterns  of  the  eigenvalues  calculated 
with  increasing  number  of  finite  elements  is  studied  carefully  and  an 
extrapolation  technique,  based  on  Richardson's  hm  deferred  approach 
to  the  limit  [8]  is  proposed. 

Recently,  various  numerical  and  analytical  solutions  related  to 
this  work  have  been  published.  Benthem  [9]  and  Kawai,  Fagitaui,  and 
Kamagai  [10]  obtained  different  analytical  solutions  to  the  problem 
of  a  Mode  I  crack  whose  front  edge  and  plane  are  perpendicular  to 
the  surface.  However,  there  is  some  question  on  the  convergence  of 
the  method  presented  by  the  latter  authors,  whose  solutions  disagree 
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with  the  results  of  Benthem  as  well  as  the  present  work.  Also, 
significant  progress  has  been  made  in  potential  theory  problems  by 
Morrison  and  Lewis  [11]  and  by  Keer  and  Farihar  [12] »  [13].  The  former 
authors  succeeded  in  obtaining  a  differential  equation  with  the  use 
of  conical  coordinates  suited  for  their  problems.  Keer  and  Parihar 
obtained  a  singular  intergral  equation  which  is  solved  numerically. 
They  extended  their  solution  to  some  three-dimensional  singu lari ties 
in  the  interior  of  an  elastic  space  which  is  irreducible  to  potential 
theory  [14a].  Keer  and  Parihar  also  solved  the  problem  of  rigid  corner 
stamp  of  small  angle  on  a  semi-infinite  body  for  which  the  solution 
is  complex  [14b] . 

These  solutions  provide  a  valuable  check  on  the  present  method 
of  solution  and  are  compared  separately  in  Chapter  III,  where  further 
numerical  results  are  presented.  These  are;  a  crack  whose  front  edge 
and  plane  are  perpendicular  to  the  surface,  where  it  is  shown  that 
Modes  II  and  III  are  coupled  and  inseparable  at  the  surface  point; 
a  crack  whose  front  edge  is  inclined  but  whose  plane  remains  perpen¬ 
dicular  to  the  surface  in  all  modes;  and  a  crack  whose  front  edge  and 
plane  are  inclined  in  all  modes.  From  energy  flux  considerations, 
these  results  show  that  upon  propagation  the  surface  point  of  a  crack 
in  symmetric  opening  will  trail  behind  its  interior,  while  in  anti¬ 
symmetric  openings  the  surface  point  will  move  ahead  of  the  crack 
interior.  It  is  also  shown  that  for  an  orthogonal  crack  in  combined 
openings,  the  surface  point  will  propagate  in  either  symmetric  mode 
or  antisymmetric  modes,  but  not  in  a  combination.  The  numerical 
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results  are  compared  with  some  recently  reported  fractographlc 
measurements  provided  by  Bell  and  Feeney  [15]. 

Recently,  John  P.  Benthem  of  Delft,  Netherlands  privately 
communicated  (April,  1978)  results  based  on  a  finite  difference 
method  directly  applied  to  the  differential  equations  of  equilib¬ 
rium  and  boundary  conditions.  His  results,  which  have  not  yet 
been  published,  are  also  compared  in  Chapter  111.  They  agree 
quite  closely  with  the  present  solutions. 


CHAPTER  I 


VARIATIONAL  EQUATION  FOR 
THE  EIGENSTATES 


1. 1  Introduction 

The  most  powerful  method  used  to  determine  near-singularity  fields 
is  that  of  asymptotic  analysis  (separation  of  variables).  The  method 
was  first  used  by  Knein  [1]  (who  thanked  T.  von  KArmdn  for  suggesting 
the  basic  approach)  in  a  problem  of  plane  elasticity,  later 
solved  independently  by  Williams  [2]  and  Karp  and  Karal  [16].  Recently, 
various  authors  have  extended  this  procedure  to  three  dimensions  to 
investigate  the  near-singularity  behavior  for  different  problems 
[3,11-14}. 

The  present  work  makes  use  of  the  same  method  in  order  to  formu¬ 
late  a  general  variational  equation  for  cracked  or  notched  linear 
elastic  bodies.  The  formulation  consists  of  the  classical  diffe¬ 
rential  equations  of  equilibrium  of  linear  elasticity  and  the  boundary 
conditions  associated  with  the  problem.  A  variational  statement  is  pro¬ 
posed  and  then  reduced  to  a  variational  equation  applicable  to  nu¬ 
merical  methods. 

The  variational  equation  determines  the  behavior  of  the  material 
near  the  point  of  singularity.  The  basic  assumption  is  that  near 
such  a  point  the  leading  terms  of  displacement  components  are 
of  the  form  r^F(0,<j>),  where  X  is  an  unknown  constant  and  F(.0,<|>)  an 
arbitrary  function  of  the  angles  0  and  <f> . 
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1.2  Governing  Equations  of  Elasticity 

The  problems  considered  in  this  work  deal  with  singularities  in 
a  linear  elastic  material.  As  a  first  step,  the  equations  governing 
such  material  are  introduced.  For  reasons  which  will  become  clear  in 
later  sections,  these  equations  are  written  in  the  spherical  coordinate 
system  (r,9,0)  and  in  terms  of  their  respective  displacement  components 
(u,v,w).  Ho  dynamic  terms  or  body  forces  are  present  since  the 
solution  is  confined  to  small  neighborhoods  of  the  singularity. 

A. )  Equations  of  equilibrium. 

The  well  known  classical  differential  equations  of  equilibrium 
expressed  in  terms  of  dilatation  and  rotation  and  transformed  to  spherical 
coordinates,  take  the  form  [Ref.  17,  pages  141  and  56]: 

—  xa  -s  3ttU 

(X  +2u)r  sin  9  H  -  2|j{:^  (ui^sin  6)  -  }  -  0  (1.1a) 

-  Suu_  ..  __ 

(X  +2n)sin  9  ff  ^  (™0  sin  9)3  *  0  (1-lb) 

_  .  3cu 

^  +  2 ^  sin  9  30  *  2u^3r  ^9^  *  “30  ■*  “  0  (1-  1c) 

where  X  and  |i  are  the  Lame' s  constants,  L  is  the  cubical  dilatation  and 
U)r,  UJq,  and  are  the  components  of  rotation: 


L  -  - -  (r2u  sin  9)  +  ^  (rv  sin  9)  +  ^  (rw)}  (1.2a) 


r2sin  9  3r 
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(1.2b)  ' 

(1.2c) 

(1.2c) 

B. )  Strain-displacement  relations. 

When  a  body  is  slightly  deformed  the  strain-displacement  re¬ 
lationships  written  in  spherical  coordinates  take  the  form  [Ref.  17, 
page  56]: 


2*r  *  { & (rw  sin  6)  ■  * (rv)} 

*  T7T  _  *  (tw  sin  9)} 

r  sin  8 

-  i  &  (»>  -  & 


e  =  u. 
rr  r 


e09  “  r  v9  +  r  U 


e00  *  nkr  w0  +  rvcot0+ru 


eftj  »  —  W_  -  —  W  COt  0  +  - 7 - rr  V 

00  r  9  r  rsin90 


1  .  1 

e,  - - — r  u  +  w  -  -  w 

0r  r  sin  9  0  r  r 


1  .  1 

e  .  -  v - v  +  —  uQ 

r9  r  r  r  9 


(l-3a) 

(1.3b) 

(1.3c) 

(1.3d) 

(l*3e) 


(l-3f) 


where,  single  subscripts  on  u,  v,  and  w  indicate  partial  derivatives, 
and  double  subscripts  on  e  indicate  strains. 

C. )  Stress-strain  relations 

When  a  linear-elastic  isotropic  homogeneous  material  is  slightly 
strained  the  stress  components  are  linear  functions  of  the  strain 
components.  With  the  strain  components  defined  by  Eqs.  (1.3  a-f)  the 
stress  components  are  [Ref.  17 f  page  126]: 


a„  -  X  e  +u.u  err 


a09  *  X  e  +Wu  e86 


°00  *  X  6  +Uu  *00 


a90  *  ^  e00 


a,  =  u  e, 
0r  M  0r 


ar0  *  u  er0 


(l.ha) 

(lAb) 

(l.hc) 

( l-4d) 

(l.he) 

(lAf) 


where. 


e 


err  +  e00 


+  e 


00 


(l.lt-s) 
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Then,  the  surface  tractions,  T^,  on  a  surface  of  unit  normal,  n^, 
written  in  index  notation,  are  [Ref.  18,  page  64]: 

Ti  "  °ijSj  *  m  (1.5) 

0.)  Comment 

It  is  well  known  that  the  material  behavior  near  a  crack  tip  or 
a  notch  apex  will  be  nonlinear.  As  is  also  typical  of  all  cases 
possessing  singularity  regions.  The  stress  values  will,  in  fact, 
become  unbounded  at  the  singularity,  although  actually  the  maximum 
stress  cannot  exceed  that  at  which  plastic  flow  takes  place.  Neverthe¬ 
less,  the  theory  of  linear  elasticity  can  adequately  describe  the 
not  too  close  stress  field  if  the  plastic  region  is  small.  The  present 
work  is  limited  to  such  cases. 

1.3  The  Williams'  Method 

Consider  the  mathematical  representation  of  a  crack  plane  inter¬ 
secting  a  semi-infinite  elastic  body.  Fig.  1.1  shows,  as  an  example, 
a  crack  plane  as  well  as  its  crack  front  edge  to  be  normal  to  the 
surface.  For  illustration  purposes,  an  imaginary  body  is  cut  out  by 
a  spherical  surface  of  small  radius  and  centered  at  the  point  where 
the  crack  front  edge  and  the  surface  meet,  i.e.,  point  0,  Fig.  1.1. 

This  normal  presentation  is  used  to  indicate  spherical  coordinates, 
although  the  equations  of  elasticity  apply  to  the  entire  body. 


A.)  Separation  of  variables. 

Consider  the  point  0,  located  at  the  smooth  line  00'  terminating 
at  point  0,  e.g. ,  Fig.  1.1.  Let  r,0,0  be  a  spherical  coordinate  sys¬ 
tem  centered  at  point  0,  such  that  ray  0=0  coincides  with  the  crack 
front  edge,  line  00*.  It  will  be  assumed  that  in  the  vicinity  of 
point  0  the  displacement  components  in  the  r,9,0  directions  are  func¬ 
tions  whose  dependence  on  r  can  be  separated  from  their  dependence  on 
9  and  0,  i.e.,  the  separation  of  variables  technique  will  be  applied: 


u(r,9,0)  «  rXF(0  ,0) 

(l.6a) 

v(r,0,0)  =  rXG(0,0) 

(1.6b) 

w(r,0,0)  =  rXH(0,0) 

(1.6c) 

with  the  restriction  that  the  exponent  1  have  a  limiting  value 
Re(X)  >  in  order  for  the  strain  energy  to  remain  finite  near 
point  0,  r  -*  0.  Hence,  the  objective  of  this  work  is  to  find  the 
smallest  possible  value  for  the  exponent  X,  Re(\)  >  giving  the 
gravest  state  of  stress  for  the  vicinity  of  point  0. 

In  this  sense,  point  0  is  considered  a  singular  point,  line  00' 
a  singular  line,  and  X  the  eigenvalue. 

The  proof  for  the  well-established  theorems  of  uniqueness  and 
existence  for  the  problems  considered  here  is  beyond  the  scope  of 


this  work.  Even  for  regular  finite  regions,  these  theorems  are  al¬ 
ready  "not  distinguished  by  simplicity"  [19,  page  89]*  Benthem  [9] 
who  solved  analytically  one  of  the  problems  presented  in  Chapter  III, 
writes:  "The  following  theorems,  though  without  proof,  will  be  con¬ 
sidered  to  be  valid  for  elastic  regions  in  the  form  of  infinite  cones 
(semi-infinite  bodies). 

i. )  If  an  infinite  conical  region  is  loaded  by  stresses  which 
behave  along  the  generators  like  r^~*  and  the  displacements  are  pre¬ 
scribed  which  are  zero  or  behave  like  r^,  then  there  is  generally  a 
solution  for  the  interior  stresses  of  the  form  axx  *  r^-1fxy(90) ,  etc. 
with  the  exception  of  an  infinite  enumerable  set  of  values  for  X. 

ii. )  For  every  value  of  X  of  the  infinite  enumerable  set  meant  under 
(!•).  there  exists  a  state  of  stress  given  by  the  above  expressions, 
whereof  the  prescribed  stresses  and  displacements  are  zero.  Such 
states  of  stress  are  called  the  eigenfunctions  of  the  cone  in  question, 

iii. )  The  infinite  enumerable  states  of  stress  (with  Re(X)  >  -§•) 
meant  in  (ii.),  are  able,  in  principle,  to  meet  every  set  of  three 
boundary  conditions  at  r  *  constant  (a  finite  cone)  provided  the 
boundary  conditions  do  not  require  a  concentrated  force  or  moment  at 
the  vertex. 

These  three  theorems  are,  if  not  proved,  generally  accepted  in 
the  corresponding  two-dimensional  analysis  of  wedges  [2,20,21]." 

B-)  Modified  equations  of  equilibrium. 

Substituting  the  expressions  (1.6a-c)  into  the  differential 
equations  of  equilibrium  (1. la-c)  it  was  found  that  the  radial  co- 
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ordinate,  mainly,  r  "  factors  out  of  the  equations.  The  following 
equations  of  equilibrium  in  the  r,9,  and  0  directions  in  terms  of  the 
functions  F,G,  and  H  and  the  exponent  X  result: 


Xr  -  (Q+2)(X  -  1)(XF  +  F  +  G0+G  cot  8  H0)  -  [(X  +  l)G0  -F00] 


-  cot  0[(X  +  1)G  -  F0]  +  .  \  (  .1  --  F  -  H  -  XH  )  =  0 
0  sin  9  sin  0  00  0  0' 


(l*7a) 


xe  ’  <G+2><XFe  +2Fe  +Gee  +Ge  coc  9  '  7^  G+nTe  H9 « 


-■M-  H.)  -  a  (Hft  .  +  H.  cot  0  -  G.  ) 

sin0  ^  Sin  ®  o0  0  sin  0  00' 


+  X[(X  +  1)G  -  F0]  =  0 


(1.7b) 


}  T  (Q  +  2)  (XF  +2F .  +  G,  cot  0  4-  — H  ) 
0  sin  0  0  0  0  sin  9  00 


VH'XB) +  (SE5+H9  cote*~^TE 

*  sm  9 


cos  9  G  .  . G  n  .  o 
sl»29  «  si”9  V 


(1.7c) 


where,  subscripts  of  F,G,  and  H  denote  partial  derivatives,  e.g. , 
2  2 

F00  ■  5  F/30  ,  v  =  Poisson  ratio, 


11* 


Q  =  2v/(l-2v) 


(1.8) 


and  X  ,XQ  »  X  symbolically  represent  the  new  modified  equations  of 
equilibrium  in  the  r,  0,  and  0  directions,  respectively. 

C- )  Modified  stress-strain  relations. 

Substituting  Eqs.  (1.6a-c)  into  the  expressions  for  the  spherical 
stress  components,  Eqs.  (1.4a-f),  the  following  modified  stress  ex- 
presssions  result: 


rr 


art  -  QftF  +  2F  +  G0  +  a  cot  0  +  ^  H0)  +  2XF 


( 1- 9a) 


Sr6  =  urX-l  ar9  "  XG  '  G  +  F0 


(l-9b) 


*  “FT  ct86  “  Q(XF  +2F  +  Ge+Gcot0+  H0) 


98 


+  2(Ge  +  F) 


se <t>  "  urx-i  a80  "  ^ 


-  H  cot  8  + 


sin  8 


(i-9c) 

( 1-9*!) 


Sr0  "  ur,X-l  cr0  =  sin  8  F0  +  XH  H 


(l-9e) 
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V  -jhr  °o» ' Q(XF  +  2F  +  S +  G  cot  6  +  nr?  V 
+  2<ilH  »4  +  G  cot  6  +  F>  U'M) 

in  which  (j  is  the  elastic  shear  modulus,  and  srrJ’‘*sS00  symbolically 
represent  the  stresses. 

D. )  Comment. 

It  is  interesting  to  note  that  one  may  expect  an  infinite 
enumerable  number  of  real  and  complex  eigenvalues  i.e.,  not  a  con¬ 
tinuous  spectrum,  for  the  problem  of  a  serai-infinite  homogeneous  body. 
Furthermore,  the  real  part  of  each  of  the  complex  roots  with  positive 
real  part  is  always  greater  than  the  smallest  positive  real  root, 
which  is  also  the  case  for  plane  problems  and  is  rigorously  proved  in 
Ref.  [16].  Therefore,  the  dominant  term  that  governs  the  behavior 
near  point  0,  Fig.  1. 1,  is  given  by  the  smallest  real  root. 

Of  course,  this  last  observation  will  not  apply  to  problems  of 
a  crack  intersecting  two-material  interfaces.  In  such  cases  all 
eigenvalues  \  are  expected  to  be  complex  in  nature,  such  as  those  found 
in  plane  strain  elasticity  [22, 2J].  The  term  of  interest  will  then 
be  given  by  the  smallest  Re(\). 

It  may  be  desired  to  write  a  generalized  Fourier  analysis  to 
include  the  displacement  fields  for  different  eigenvalues  which  must 
result  from  the  solution  of  an  infinite  set  of  equations.  However,  the 
orthogonality  properties  in  three-dimensions  might  be  insufficient  to 
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determine  their  participation.  This  is  also  the  case  in  plane  problems, 
as  Williams  [24]  noted:  "...  let  it  suffice  to  point  out  that  the 
solutions  f(\,a,v)  *  0  yield  an  infinite  number  of  eigenvalues  which 
may  be  complex.  After  ordering  these  values  according  to  their  absolute 
magnitudes,  one  may  construct  from  them  an  infinite  set  of  eigenfunc¬ 
tions  whose  elements  are  non-orthogonal.  Furthermore,  the  completeness 
of  the  set,  although  intuitively  probable,  has  not  been  mathematically 
established."  The  present  literature  seems  to  lack  the  proof  for  the 
completeness  of  the  set,  even  for  two  dimensional  problems. 

In  any  event,  let  such  set  intuitively  exist  for  the  three 
dimensional  problem.  Then,  in  the  vicinity  of  point  0,  the  behavior 
of  the  displacement  and  stress  fields  will  be  determined  by  the  eigen¬ 
function  characterized  by  the  smallest  eigenvalue;  provided  the  loading 
near  such  point  is  not  critical,  see  theorem  (iii.)  page  12. 

1.4  Construction  of  a  Variational  Equation 

A.)  Cartesian  (9-<j) -plane. 

Finite  element  studies  are  simplified  when  flat  planes,  rather  than 
curved  surfaces  are  used.  Expecting  that  the  solution  to  these  problems 
will  make  use  of  numerical  techniques,  the  domain  O'ACO  of  a  unit 
sphere  from  Fig.  (1.1)  is  visualized  in  a  fictitious  (0-0)-plane  shown 
in  Fig.  (1.2).  This  approach  has  been  successfully  developed  in 
Ref.  [3, page  226]. 

The  singularity  ray  00',  Fig.  (1. 1)  placed  on  the  pole  of  its 
spherical  coordinate  system  appears  in  the  (0-0) -plane  as  a  straight 
line  segment  at  0  »  0.  The  surface  of  the  semi-infinite  body  6  ■  tt/2, 
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0  £  0  £  2tt  appears  in  the  (0-0) -plane  as  a  straight  line  segment 
0  *  tt/2,  0*0*  2tt,  etc. 

Let  then  n  =  (nQ,n  )  be  the  unit  normal  to  the  surface  of  the 
—  0  0 

body  when  plotted  in  the  (0-0) -plane,  with  0  and  0  being  regarded  as 

the  cartesian  coordinates  in  such  a  fictitious  plane.  Thus,  n  = 

(-d0/ds,d0/ds)  where  s  is  the  length  of  a  boundary  curve,  or  n  /n  = 

9  0 

-d@/d0  where  d0  and  d0  are  increments  along  such  boundary  in  the  0 
and  0  directions,  respectively. 


B.)  Free  surface  conditions. 

It  will  be  assumed  that  in  a  sufficiently  small  neighborhood  of 
point  0,  Fig.  (1.1),  there  are  no  loads  applied  at  the  body  surfaces 
(surfaces  formed  by  radial  rays  emanating  from  point  0) ,  or  at  the 
cracked  surfaces. 

Indeed,  the  purpose  of  this  work  is  to  determine  all  possible 
states,  called  eigenstates  which  are  strictly  characterized  by  the 
eigenvalues  Hence,  surface  loads  may  be  prescribed  at  body  surfaces 
sufficiently  remote  from  point  0,  for  which  the  differential  equations 
of  equilibrium  but  not  the  boundary  conditions  will  be  satisfied. 
According  to  the  principle  of  . superposition  the  actual  state  of  stress 
for  given  boundary  conditions  can  be  expressed  as  a  linear  combination 
of  its  eigenstates.  But,  as  far  as  the  eigenvalues  X  are  concerned 
these  boundary  conditions  are  irrelevant,  as  it  is  also  true  in  planar 
problems,  where  X  =  \  for  any  loading  combination. 

Therefore,  without  loss  of * generality,  it  will  be  assumed  that  all 
body  surfaces  are  formed  by  radial  rays,  and  that  the  boundary  conditions 
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at  the  free  surfaces  in  close  vicinity  of  point  0,  Fig.  (1-1),  are: 


body  surface:  *  Cgr  =  *  0  at  0  *  tt/2  (1.10a) 

crack  plane:  a ^  -  a^r  =  *  0  at  0  *  0,2rr  (1.10b) 


these  boundary  conditions  written  as  surface  tractions  p  ,p  ,  and 

r  y 

p^  in  the  (8-0) -plane  take  the  form: 


Pr  *  Sr0ni0  sin  9  +  Sr0° 0  =  0 


p8  "  s09nB  Sin  0  +  S00n0  “  ° 


p .  *  s  n  sin  0  •+•  s  ..n  «  0 

0  00  0  00  0 


(1.11a) 

(1.11b) 

(l-llc) 


No  boundary  conditions  are  specified  at  infinity,  r  -*  <*;  only  the 
local  problem  of  stress  singularity  at  point  0,  r  -*  0,  is  considered, 
and  hence,  a  small  domain  about  point  0  is  required. 

C. )  Variational  statement. 

The  differential  equations  (l-7a-c)  together  with  the  boundary 
conditions  (1. lla-c)  may  be  combined  to  form  the  following  variational 
statement  in  the  (8-0) -plane: 

JJ  (Xr6F  +  XgdG  +  X$5H)sin  0d0d0  -  j  (pr5F  +  p05G  +  p^5H)ds  =  0 


(1.12) 
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in  which  s  is  the  length  of  the  boundary  of  the  region  of  the  (0,0)- 
plane;  A  is  the  area  of  such  region;  and  <5F,6G,  and  6H  are  arbitrary 
continuous  functions  of  0  and  0  which  have  piece-wise  continuous 
derivatives  and  satisfy  all  displacement  boundary  conditions,  if  any. 
Conversely,  from  the  fact  that  Eq.  (1.2)  must  hold  for  any  kinematically 
admissible  functions  6F,dG,  and  <5H  it  follows  that  Eqs.  (1.7a-c)  and 
(1. lla-c)  must  be  satisfied.  Thus,  the  variational  statement  (1.12) 
must  be  equivalent  to  Eqs.  (1.7a-c)  and  (1. lla-c). 

The  variational  statement  presented  in  Eq.  (1.12),  is  analogous 
to  the  three-dimensional  statement  one  forms  in  order  to  obtain  the 
strain  energy  [18],  except  that  integration  with  respect  to  r  has 
already  been  carried  out  in  the  unit  sphere,  since  the  r  dependence  can 
be  factored  out,  see  Eqs.  (1.7a-c)  and  (1. lla-c). 

D. )  Variational  equation. 

The  variational  statement,  Eq.  (1.12)  involves  second  derivates 
of  F,G,  and  H,  which  are  contained  in  the  expressions  for  Xr,Xg,  and 
X^,  Eqs.  (1.7a-c).  Since  numerical  techniques  give  rise  to  larger 
error  for  higher  order  derivatives,  it  is  necessary  to  transform 
Eq.  (1.12)  to  a  form  which  involves  no  higher  than  first  order  de¬ 
rivatives  of  F,G,H  and  of  5F,dG,6H.  Also,  to  be  able  to  apply  the 
finite  element  method  it  is  necessary  that  during  this  transformation 
the  boundary  integral  in  Eq.  (1.12)  be  included  in  the  variational 
equation;  otherwise  the  natural  boundary  conditions  would  not  be 
satisfied  when  the  finite  element  method  is  used. 


Indeed,  a  transformation  by  Green's  integral  theorem  [18,  page  279] 
applied  in  the  Cartesian  (6-0)-plane  has  been  found,  such  that  both 
objectives  are  reached  simultaneously.  The  formulation  is  given  in 
detail  in  Appendix  A.  The  resulting  variational  equation  is: 


JI(VF  +  +  Vr#  +  +  V°e  +  *5  %  +  V 

A  0  0  0  0  ^ 


+  ^  fiH^)sin  0  d0d0  =  0 


(1-13) 


in  which  SFg  *  3<$F/3<S, ...»  =  36H/3<j>  and  the  following  notations  are 


4p  *  CQ(  1  -  X)  +  2][(\  +2)F  +  G0  +  G  cot  0  +  H^j 


-  2X(\+2); 


$  =  (X  -  1)G  +  Fq;  A 


h  K  +  a-DHl; 


8’  *F,  sin  8  *-sin  9  <*> 

0 


%  -  {(Q+2)[(X+2)F  +  Gfl+G  cot  9+  j~  H  ]  -  2(Gq  +F) 


2XF}cot  0  -  2(F0  -G)  -X(X  +  1)G-XF0; 


%  -  Q[(\+2)F  +  Gg+G  cot  H,]  +2(Qg  +T); 
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%  -  dr?  (He  -  H  cot  6  +  dr  V 1 


+  X»+1>H+  d?  V* 


4>  »  Ho  -  H  cot  0  +  Q  G . ; 

Hq  «  sm  0  0 


ft 


H. 


d? tQ[(x  +2)F  tVe“t,tsf  V 


2(d?  V""1*  f,] 


(i-») 


where  <^,...,4^  are  not  partial  derivatives  of  some  function  4,  and 

0 

are  used  only  for  notation. 

Thus,  the  variational  statement  of  the  problem  is:  Functions  F, 

G,  and  H  are  the  solution  of  the  problem  if  and  only  if  they  satisfy 
Eq.  (1.13)  for  any  kinematically  admissible  variation  &F,&G,  and  6H. 

Existance  of  the  variational  equation  which  contains  no  boundary 
integral,  Eq.  (1. 13),  indicates  that  natural  boundary  conditions, 

Eqs.  (l.lla-c),  will  be  automatically  fulfilled  when  the  finite  element 


method  is  used. 
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£. )  Comment. 

Alternatively,  it  is  possible  to  derive  Eq-  (1-13)  from  Eq.  (1.12) 
by  means  of  Stokes  theorem  applied  to  a  unit  sphere.  It  has  been 
checked  that  this  gives  the  same  result.  It  has  been  also  checked 
that  Eq.  (1.13)  can  be  transformed  back  to  Eq.  (1.12)  by  means  of 
Gauss  or  Stokes  theorems.  For  the  sake  of  brevity  these  derivations 
are  not  given.  Instead,  an  alternate  and  independent  derivation  of 
the  one  just  given  is  derived  in  the  next  section. 

1*5  Alternative  Derivation  of  the  Variational  Equation 

The  basic  variational  equation,  Eq.  (1-13)  can  also  be  derived 
from  the  principle  of  strain  energy.  The  derivation  is  more  direct 
but  involves  certain  steps  which  were  difficult  to  foresee  at  the 
early  stages  of  this  project  without  recourse  to  the  derivation  just 
presented. 


A. )  Principle  of  minimum  potential  energy. 

The  total  potential  energy  stored  within  a  linear  elastic  body 
of  volume  V  and  surface  S,  when  no  body  forces  or  dynamic  terms  are 
present,  is: 

U  -  I  i|f  r2sin  6  drdSdtfi  -  f  (T u  +  T,v  +  T,w)ds  (1.15) 

d  V  *  g  r  “  V 

where  Tr,T0,T  are  the  surface  tractions  defined  in  Eq.  (1.5),  and  ty  is 
the  strain  energy  density: 
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TO  1  1  O  <■* 

ft  *  ^  f Q(u  +—  v0  +—  u  + - : — —  w  +  -  v  cot  6)  +  2[u  c 

’  2  r  r  0  r  r  sin  9  <f>  r  7  L  r 


+  (j  v6+i  u)2  +  v  cot  8  +  i  u)2] 


+  (v  i  v  +  i  Ue)2  t(I,e  .  i  „  cot  8  +  v  )‘ 


,  ,  1  1  s2, 

+  ( - - s-  U  +  W - w)  J 

'r  sin  8  4  r  r 


(1.16) 


in  which. 


0  =  2v/(l  -  2v) 

(1. 17a) 

4  =  E/2(l  +  v) 

(1.1Tb) 

E  is  the  Young's  modulus  and  v  the  Poisson's  ratio  characterizing  the 
linear  elastic  material. 

According  to  the  principle  of  minimum  potential  energy,  the  state 
of  equilibrium  is  a  state  for  which  the  first  variation  of  the  total 
potential  energy  vanishes.  Thus,  consider  the  displacement  variations 

6u  =  c  «  (1. 18a) 


Sv  ■  e  v 


(1. 18b) 


2k 


Sw  *■  c  w  (1* 18c) 

where  E  Is  a  variable  parameter  and  u,v,w  are  any  chosen  displacement 
distribution  which  are  sufficiently  smooth  and  satisfy  all  kinematic 
boundary  conditions.  Then,  the  state  of  equilibrium  is  determined 
by  the  first  variation  of  U,  Eq.  (1.15): 

&U  -  e[f“]  -  f  (1.4 u  +  A  6u  +  ...  +  t  6W  )r2sin  0  drd0d0 

d  e=0Jv  r  0  0 

-  f  (T  Su  +  T„bv  +  T  4w)ds  =  0  (1.19) 

-  r  o  0 


B- )  An  unorthodox  step. 

o\ 

If  Eqs.  (1.6a-c)  were  substituted  directly  into  Eq.  (1.19)  t 

would  factor  out.  However,  the  remaining  expression  would  not  be 

able  to  satisfy  the  equations  of  equilibrium  when  Gauss  theorem  is 

applied  to  Eq.  (1.19).  To  circumvent  this  critical  problem,  consider 

the  terms  <jf  ftu  ,  ;jfv  6v  ,  and  ,6w  separately  from  Eq.  (1.19)- 
r  r  r  r  r  r 

These  terms  can  be  simplified  by  Gauss  theorem  as: 

J  (*u  4u  +  hv  6v  +  i|f  6 w  ) dv  - 

“  V  ur  r  vr  *  r  r 

(K  6u)  '  +£  (*v 6v)  ■  *v 

V  r  r  r  r 

+  £  <*w  *w)  •  *w£*w  ^dV 
r  r 


(1.20) 


where  d/dr  is  the  total  derivative  with  respect  to  r. 

Applying  the  Gauss  theorem  to  the  sum  of  the  first,  third  and 
fifth  terms  of  Eq.  (1-20)  and  placing  them  back  into  Eq.  (1-19)  gives 

6U  "  J  t<*u  •  £  ♦«  )&u  +  VSue  +  V5u0 

V  3T  0  (p 

+  (*v  '  £  ^v  )iv  +  V6ve  +  V6v0 

r  0  0 

+  (*w  -  £  *w  )6w  +  Viw9  +  *w  6w0}dV 
r  0  0  r 

+  J  £(+„&«  +  %  +  *w  fw)nr 

J  S  r  r  r 

-  Tr&u  +  T06.v+T06w}ds  =0  (1.21) 

For  the  particular  problems  studied  here,  the  surface  integral  of 
Eq.  (1.21)  vanishes  because  of  the  boundary  conditions  conveyed  by 
these  problems,  i.e. : 

If  the  surface  tractions,  Eq.  (1-5) ,  are  expanded,  one  obtains 
that  on  the  free  surfaces  near  point  0  [see  Fig.  (1.1)  and  Eqs.  (l.lOa-b)J 

i)  ^  ,  which  is  not  present  on  any  surface  formed  by 

rays  emanating  from  point  0. 

it)  ar0  •  i(iv  m  0  on  the  body  surface,  0  =  tt/2.  ,  cr0  is  not  pre¬ 
send  on  the  crack  plane  0  •  0. 
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iii)  a  .  *  ijf  *  0  on  the  crack  plane  0  ■  0,  but  not  present 

r0  wr 

on  the  body  surface  0  *  tt/2. 

Iv)  because  of  the  free  surface  conditions  the  remaining  terms 
are  also  zero  and  because  of  (i-iii): 


Tr  =  °r0n9  +  V*  "  ° 


T0  “  CT00n0  +  CT00Q0  *  0 


T0  =  %Q0  +  C00n0  =  0 


(1.22a) 

(1.22b) 

(1.22c) 


which  are  analogous  to  the  boundary  conditions  given  by  Eq-  (l.lla-c). 

If  the  expressions  (1.6a-c)  are  now  substituted  into  Eq.  (1.21) 
one  would  obtain  the  same  basic  variational  equation,  Eq.  (1.13), 
after  integration  with  respect  to  r  is  performed  on  the  unit  sphere 
and  transformed  to  the  (9-0) -plane. 

C. )  Lack  of  symmetry  and  non-existence  of  a  minimum  principle. 

It  is  particularly  noteworthy  that  the  integrand  of  Eq.  (1.13)  or 
Eq.  (1.21)  is  non-symraetric,  and  so  is  the  system  of  linear  equations, 
Eq.  (1.12),  which  is  Eq.  (1.13)  applied  to  finite  elements,  (k^  ^  k^). 
This  means  that  the  variational  equation  cannot  be  written  in  the  form 
of  a  classical  stationary  principle  [18],  6W  -  0  (or  minimum  principle, 
W  *  min.),  which  would  yield  Eq.  (1.15).  For  an  elastic  material  this 
might  seem  surprising.  However,  a  deeper  analysis  indicates  that  it 


must  be  so. 
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Assume  that  the  integrand  of  Eq.  (1.13)  is  symmetric  with  F,G,  and 
H.  Then  the  discrete  eigenvalue  problem  for  X  resulting  from  Eq.  (1. 13) 
would  have  to  be  a  symmetric  matrix  *  k^,  Eq.  (2.12).  This 
implies  that  all  roots  X  would  have  to  be  real.  But  this  cannot  be 
possible  because  the  same  variational  equation,  Eq.  (1.13),  must  hold 
also  for  plane  strain  problems  with  two  material  interfaces,  whose 
solution  are  known  to  exhibit  oscillating  singularities  [22,23]  for 
which  X  is  complex.  Hence,  Eq.  (1.13)  cannot  be  symmetric.  This 
contrasts  with  the  analogous  potential  theory  problem  for  which  a 
minimum  variational  principle  in  the  (9-0) -plane  does  exist  [3 ],  with 
the  consequence  that  in  potential  theory  the  eigenvalue  X  is  always 
real. 

To  prove  that  the  variational  equation  Eq.  (1.13)  must  be  non- 
symmetric,  it  is  sufficient  to  show  that  it  must  be  so  in  the  special 
case  of  plane  elasticity.  This  can  be  done  by  dropping  out  the 
integration  over  9  and  substituting  9  *  rr/2  and  then  setting  G  *  v  =  0 
in  Eq.  (1.13).  In  that  case,  the  most  general  quadratic  functional 
involving  F(0) ,H(0) ,F' (0)  *  3F /;>0,  and  H'(0)  =  3H/&0  is: 

W  *  J^1[|(A1F2  +  A2f'  +  A^H2  +  A^H'2)  +  A^FF'  +  AgFH 

+  A^FH'  +  AgF'H  +  A^F'H*  +AloHH']d0  (1-23) 

The  associated  Euler  equations  [13]  are: 
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AjF  -  ^F”  +  AgH  +  (A7-Ag)H'  -  A9H”  «  0 

(1.23a) 

A^H  -  A^H"  +  F  -  (At-A8)F'  -  A9F"  =  0 

(1.23b) 

and  the  corresponding  natural  boundary  conditions 

are 

at  0  =  0  or  0  ~0, 

A^F'  +  A$F  +  AgH  +  A^H'  =  0 

(l.24a) 

A^H'  +  A^F  +  A^F'  +  A1QH  =  0 

(1.24b) 

The  actual  differential  equations  for  F  and  H,  as  obtained  by  sub¬ 
stituting  u  *  r^F(0)  and  w  =  r^"H($)  into  the  planar  differential 
equations  of  equilibrium  in  the  polar  coordinate  system  (r,0),  given 

by  Karp  and  Karal  [16],  have  the  form 

C^a  F  +  a2F"  +  a1H’)  =  0 

(1.25a) 

C2(boH  +  bgH"  +  bjF')  -  0 

(1.25b) 

and  the  actual  boundary  condl  ions 

CjCcjF'  +  c2H)  =  0 

(1.26a) 

W'  +  c^F)  =  0 

(1.26b) 

2$ 


where  CpCgjC^,  and  are  aribtrary  non-zero  constants;  and 
ao,a^,a2;  bo,b^,b2;  C1,C2,C3,C4  are  certa*n  given  constants. 

Equating  the  coefficients  of  all  corresponding  terms  of  Eqs.  (1.23a-b) 
and  (1.24a-b)  with  Eqs.  (1.25a-b)  and  (1.26a-b)  one  obtains  a  system  of 
14  linear  algebraic  equations  for  Ap...,AjQ,  Unknowns 

A^,...,A^q  can  easily  be  found,  which  leaves  a  system  of  four  linear 
equations  for  C^,...,C^  which  are  homogeneous.  The  determinant  of  this 
equation  system  was  found  to  equal  Because  X  cannot  be  restricted 
to  equal  zero,  it  follows  that  C^,...,C^  cannot  be  non-zero.  Thus, 
there  is  no  way  to  make  Eqs.  (1.23a-b)  and  (1.24a-b)  equivalent  to 
Eqs.  (1.25a-b)  and  (1.26a-b),  which  means  that  a  variational  functional 
W  does  not  exist  for  plane  problems.  So,  it  cannot  exist  for  the 
three-dimensional  problem  as  well. 

D.)  Comment. 

Indeed,  the  two  derivations  take  somewhat  different  procedures 
from  those  found  from  classical  variational  methods  in  linear  elasticity. 

But,  it  is  noteworthy  that  both  derivations  complement  each  other  in 
the  following  manner: 

a. )  The  boundary  conditions  in  the  (9-0)-plane,  Eqs.  (l.lla-c) 
have  to  be  included  in  the  first  derivation  (see  Appendix  A),  while 
in  the  second  they  vanish  because  of  the  free  stress  surfaces,  Eqs. 

(1.22a-c). 

b. )  In  the  first  derivation  the  boundary  conditions  (mainly,  c? 

not  being  present  on  any  surface  when  surface  tractions  are  considered;  i 

I 

■1 
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<jrg  *  0  on  9  *  rr/2,  but  not  present  on  0  ■  0;  ar(^  =  0  on  0  *  0,  but  not 
present  on  9  *it/2),  allows  the  formation  of  the  (9-0)-plane,  other¬ 
wise,  the  radial  component  of  the  unit  normal  vector,  n^,  on  the 
surfaces  would  be  present  in  Eq.  ( 1. lla-c).  In  the  second  derivation 
they  come  out  as  a  result  of  the  application  of  Gauss  theorem,  but 
vanish  due  to  the  reasons  just  presented. 

c.)  The  first  derivation  makes  use  of  the  differential  equations, 
while  the  second  starts  from  the  total  potential  energy.  Both  pro¬ 
cedures  end  with  the  same  variational  equation  after  some  intuitive 
manipulations. 

The  lack  of  symmetry  and  the  non-existance  of  a  minimum  principle 
for  the  variational  equation  corresponds  to  the  fact  that  both  real 

and  complex  eigenvalues  must  be  associated  with  this  variational 
equation. 

It  must  be  stressed,  however,  that  the  entire  present  formulation 
is  contingent  upon  the  assumption  of  the  separated  form  of  the  eigen¬ 
state  (Eq.  1.6,  page  10).  There  exists  no  proof  that  the  eigenstate 
ought  to  have  this  form,  and  that  other  eigenstates,  possibly  even  not 
separated  ones,  might  exist  and  might  be  more  severe,  even  though  this 


seems  unlikely. 


CHAPTER  II 


METHOD  OF  SOLUTION  ON  A 
FINITE  ELEMENT  GRID 


2. 1  Introduction 

The  variational  equation  derived  in  Chapter  I  and  given  by  Eq.  (1.13) 
has  the  tremendeous  advantage  that  the  stress  boundary  conditions  are 
automatically  implied  whenever  a  free  surface  is  present.  Therefore, 
compared  to  a  finite  difference  method  where  free  surfaces  would  require 
additional  programming,  the  finite  element  method  is  selected  to  ap¬ 
proach  the  problem.  The  program  is  then  written  in  the  Fortran  IV 
computer  language. 

The  finite  elements  are  chosen  as  simple  four-node  quadrilaterals. 

The  distribution  functions  for  F,G,  and  H,  Eq.  (1.6a-c),  are  considered 
bilinear  in  9  and  0.  The  coefficients  of  the  stiffness  amatrix, 

Eq.  (2.9)5  ate  calculated  by  the  Gaussian  numerical  integration  tech¬ 
nique  using  nine  integration  points,  [25]. 

The  variational  equation  emerges  as  a  generalized  non-linear  problem 
for  the  eigenvalues.  Various  methods  of  numerical  solution  of  this 
type  of  problems  have  been  discussed  in  detail  in  Ref.  [3].  Method 
B  from  page  230  of  Ref.  [3]  has  been  selected  to  search  for  the  root  X. 

A  method  of  solution  when  X  is  complex  has  also  been  discussed  in 
Ref.  [6]  in  connection  with  other  problems.  The  root  of  smallest 
value,  or  of  smallest  Re(X)  in  the  case  of  complex  roots,  is  of  main 
practical  interest. 
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An  extrapolation  technique  based  on  the  "deferred  approach  to 
the  limit"  [8]  is  proposed  for  the  final  value  of  root  X  as  the 
number  of  grid  subdivisions  goes  to  infinity. 

The  program  is  general  and  capable  of  handling  various  situations, 
such  as  intersections  of  crack  plane  and  crack  front  edge  of  any 
orientation,  notches  of  any  orientation  and  of  any  opening,  etc.  The 
program  will  be  also  capable  of  handling  cases  when  the  exponent  X  is 
expected  to  be  complex,  such  as  intersections  of  crack  edges  with  two- 
material  interfaces.  This  requires  a  conversion  of  the  Fortran  pro¬ 
gram  to  complex  arithmetic. 

2.2  Finite  Element  Formulation 

A.)  Treatment  of  line  singularities. 

From  the  three-dimensional  singularity  point  0,  Fig.  1.1,  page  11, 
there  usually  emanates  a  stress  singularity  line,  such  as  the  crack 
front  edge  shown  as  line  00'  which  coincides  with  the  polar  ray 
9*0.  The  displacements  near  this  line  usually  behave  as  (rp)p  [3,11] 
such  that  rp  represents  the  distance  from  the  ray  9=0  when  9  -*  0. 

The  exponent  p  will  then  represent  the  exponential  behavior  for  the 
displacement  field  near  the  singularity  line.  For  the  crack  front 
edge  considered  here,  the  values  p  =  0,^-, 1,3/2,. . . ,  are  possible,  and 
for  notch  edges  other  values  of  p  would  apply. 

From  the  theory  of  the  finite  element  method  for  plane  problems 
it  is  known,  for  example,  that  the  rate  of  convergence  in  the  pre¬ 
sence  of  square-root  singularity  is  only  0(h),  while  in  its  absence 

2 

the  convergence  is  quadratic,  0(h  ),  h  being  the  element  size  [26]. 
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It  is  conceivable  that  the  functions  F,  G,  and  H,  Eq.  (1.6a-c), 
may  exhibit  gradient  singularities  at  the  point  where  the  crack  front 
edge  00*  (a  singularity  line)  emanates  from  point  0,  Fig.  1.1;  then 
p  <  1.  Such  functions  are  not  suitable  for  numerical  calculations, 
and  if  they  are  approximated  numerically,  their  accuracy  and  conver¬ 
gence  are  adversely  affected  by  the  presence  of  singularities.  This 
difficulty  can  be  avoided  by  using  singular  finite  elements  near 
the  singularity  line.  A  more  convenient  method  has  been  proposed  and 
used  with  success  in  Ref.  [}].  In  this  method,  the  displacements  in 
the  r,0,  and  0  directions  are  expressed  as 


u(r,9 ,0)  »  rI1r1Pf(0,0)  =  r^pPf(9,0) 

(2.1a) 

v(r,9,0)  *  ^^(0,0)  *  rXpPg(9,0) 

(2.1b) 

w(r,9,0)  *  rnr^ph(9 ,0)  =  r^pph(9,0) 

(2.1c) 

in  which,  p  is  the  exponent  for  the  displacement  field  near  the  singu¬ 
larity  line;  X  *  n + p;  r^  =  rp;  p  is  any  chosen  smooth  continuous 
function  of  9  and  0  which  is  non-zero  everywhere  except  on  the  singu¬ 
larity  line  9  ■  0,  and  which  represents  the  distance  measured  on  a 
unit  sphere.  Possible  choices  are  p  *  0,  p  ■  sin  9,  etc.  The  second 
suggestion  will  be  used  for  numerical  calculations,  since  p  =  sin  9 
will  then  represent  the  exact  distance  from  the  ray  not  only  for 
9  -*  0,  but  everywhere  in  the  domain.  (Note,  however,  that  p  ■  sin  9 
cannot  be  used  when  the  angle  9  ■  tt  is  part  of  the  domain  and  where 


3^ 


no  line  of  singularity  exists,  i.e.,  at  6  *  rr)-  Thus,  it  is  con' 
venient  to  introduce  the  notations: 


F(9,0)  -  pP  f(9,0) 

(2.2a) 

0(0,0)  *  pp  g(e,0) 

(2.2b) 

H(0,0)  -  pP  h(9,0) 

(2.2c) 

PP  =  (sin  8)p 

(2. 2d) 

If  the  field  near  the  singularity  line  varies  as  p 2  when  p  is  set  to 
then  functions  f,g,  and  h  may  be  expected  to  be  free  of  gradient 

singularities.  This  would  make  the  convergence  rate  quadratic, 

2  10 
0(h  ),  [26].  On  the  other  hand,  if  components  of  types  p  ,p  ,  and 

possibly  other  components  of  different  exponents  were  present  in  the 

solution  [9],  the  rate  of  convergence  would  not  be  quadratic,  but 

slower  than  quadratic  [26]. 

When  several  exponents  p  are  p..  sent,  the  lowest  one  must  be  used. 
This  is  shown  as  follws.  Consider  that  an  exponent  p*,  which  differs 
from  the  actual  value  of  p,  is  used.  Then,  the  displacement  and  stress 
fields  would  behave  as: 

ut  ~  r^  9P  F(8 ,<t>) 


(2-3a) 
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o*ij  «■*  ou/30  <•»  r^8P  ^F(8 j0)  (2.3b) 

for  the  exact  solution,  and 

* 

Ui  ~  rX  9P  F*(0,0) 

* 

~  au/39  -  rXeP  "V(0,0) 

for  the  numerical  solution. 

Equating  the  two  critical  expressions,  for  the  stresses  a^, 
one  obtains : 

♦ 

F*(0,0)  *  0P’P  F(0,0)  (2.5) 

* 

where  the  function  F(9,0)  is  bounded.  If  one  chooses  p  >  p,  function 

* 

F  (0,0)  can  obviously  become  unbounded  as  0  ->  0;  and  it  cannot  be 
adequately  represented  numerically.  Thus, 

p  <  p  (2.6) 

* 

is  necessary.  The  best  choice  would  naturally  be  to  make  p  equal  to 
the  lowest  exponent  p  present. 

Under  these  considerations,  there  is  still  the  restriction  that 
along  the  crack  front  edge,  9  -»  0,  the  displacement  field  must  exhibit 


(2. ha) 

(2.hb) 


the  behavior 
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9^,  such  chat  q  >  0  (2.7) 

for  the  strain  energy  to  remain  bounded.  Note  that  q  is  not  the 
exponent  for  the  singularity  function  p. 

B. )  Displacement  distribution  in  an  element. 

Choosing  a  finite  element  grid  in  the  (8-0) -plane,  Fig.  2.1, 
the  unknown  functions  F,G,  and  H  may  be  represented  within  each  finite 
element.  Fig.  2.1,  in  the  form 

M 

F(9,0)  X.F1  ,  F1  -  pPfi(0,0)  (2.8a) 

i-1 

M 

6(0,0)  \Gl  ,  G1  -  pPgi(0,0)  (2-3b) 

i-1 

M 

H(0,0)  XiHl  ’  h1  *  PPhi(9,0)  (2-3c) 

i-1 

in  which  X i  -  1,2, are  the  nodal  values  of  f(9,<f>),  g(8,0)  and 
h(0,0);  and  f*(8,0),  gi(0,0)  and  hi(0,0)  are  given  distribution  functions 
within  the  finite  elements,  usually  chosen  as  polynomials  in  0  and  0 

[25]. 

Denoting  (8  ,0  )  as  the  coordinates  of  the  node,  the  dis¬ 
ci  m 

tribution  functions  must  be  chosen  such  that 
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(2. 10c) 


(2.H0 


kij  -  JK rJ  +<re +  *iA +Acl*  4A +  iA +  ,i'hhJ 

^  y  (p  y  (p 

+  $H  Hq  +  H0^  Sitl  6  d8d0 

6  0 

Note  that  the  stiffness  matrix  [k^.]  is  non-symmetric;  i.e.,  k^  ^  k^ 
in  general.  The  variational  equation  (2.12)  must  hold  for  any  choice 
of  6X^  (i  =  and  this  requires  that 

M 

Y  kijXj  =0  (i  -  1,...,M).  (2.15) 

i*l 

This  is  a  system  of  M  linear  homogeneous  algebraic  equations,  re¬ 
presenting  an  eigenvalue  problem.  All  stiffness  coefficients  k„,  not 
just  the  diagonal  ones,  depend  on  singularity  exponent  X,  and  so  the 
eigenvalue  problem  is  of  the  generalized  type.  Furthermore,  it  is 
easy  to  see  that  are  polynomials  in  X,  as  well  as  in  Poisson  ratio 
v  (when  multiplied  by  l-2v); 

kij  *  kij(X,v^  (2-16) 


D. )  Integration  for  the  stiffness  matrix. 

The  finite  elements  are  chosen  as  four-node  quadrilaterals  in  the 
(0-$)-plane.  Three  degrees  of  freedom  are  placed  at  each  node  in  order  to 


it-o 


accomodate  the  displacement  field  (f,g,h).  The  basic  distribution 
shape  function  fi(8,0),  g*(9,0),  and  h*(0,0)  on  the  original  rectangle 
are  considered  as  bilinear  in  8  and  0,  i.e.,  a  +  b0  +  c0  +  d00*. 

Following  the  conventional  methods  found  in  finite  element  tech¬ 
niques  [25 ,27], the  finite  element  stiffness  matrix  is  obtained  by 
mapping  a  general  quadrilateral,  Fig.  ^2..1)}  into  a  unit  square. 

Fig.  (2.2;,  given  by  the  transformation 


!  0 


,  0 


~  > 


1 


(2. 17) 


F  ’ 

'v. 


where 


B  =  (B1,B2,B3,B1+)  (2.18a) 

B^  =  i(l+8  9i)(l+0  0 j ) ,  i  =  k, l,m,n  (2.18b) 

I  =  (0K’eL’0M’eN)X  (2>l8c) 

l  ■  (0k)0l50m50n)T  (2.l8d) 


in  which  the  subscript  T  denotes  a  transpose;  (8^,0^)  are  the  corners 
(+  1,+  1)  of  the  unit  square  numbered  clockwise  beginning  at  (-1,-1); 
(8l,0j)  are  the  corresponding  corner  coordinates  of  the  quadrilateral 

_ 

In  retrospect,  it  appears  that  much  more  accurate  results  could  have 
been  obtained  with  higher-order  finite  elements. 
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*■  $ 


Fig.  2. 1:  Finite  element  grid 
in  the  (0-0) -plane. 


(-1,1) 


6* 


Fig.  2.2:  Unit  square,  obtained 
by  mapping  a  general 
element  from  Fig.  2.  1 , 
using  Eq.  (2.  IT). 


(1,1) 
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element;  (9  ,0  )  are  the  coordinates  of  a  general  point  within  the 
unit  square;  and  (9,0)  are  the  coordinates  of  the  corresponding  point 
on  the  quadrilateral  element. 

Carrying  out  the  foregoing  transformation  of  variables  0  and  0 

*  * 

into  9  and  0  ,  the  stiffness  coefficients  given  by  Eq.  (2.14)  may  be 
expressed  in  the  well-known  manner  [25]: 


JJ  Y(0,0)d0d0  =  J1  j1  Y*(9*,0*)d0* 


-1  -1 


d0 


(2. 19) 


in  which  Y(9,0)  is  the  integrand  of  Eq.  (2.14);  and  where 


J  = 


,0  )  »  |J| 

*(e  ,0  ) 

(2.20a) 

f 

.1 

'  ae/ae' 

,  a0/ae 

(2.20b) 

'  39/30^ 

,  30/30*  1 

J  being  the  Jacobian  of  the  transformation  given  in  Eq-  (2.16). 

The  integration  of  the  second  term  in  Eq.  (2.18)  is  carried  out 
numerically  by  the  Gaussian  quadrature  formula  [25],  over  nine  points 
a^  of  weights 


ij 


3 

e  * 

L 

p*l 


3 

r7 

L 

q»l 


H  H 

p  q 


(a  ,a  ) 
v  p*  q 


(2.21) 


^3 


in  which  a^  *  -a^  =  0.774596669241493;  =  0;  =  5/9; 

Hg  =  8/9- 

A  (12  x  12)  stiffness  matrix  is  thus  obtained  for  all  elements, 
which  are  then  incorporated  element  by  element  into  the  final  global 
stiffness  matrix  k„.  A  detailed  Fortran  program  is  given  in 
Appendix  B. 

E- )  Comment. 

The  treatment  of  line  singularities  described  in  Section  2.2A 
may  be  thought  of  in  a  different  light.  The  representation  of  the 
displacements  with  the  approximate  form  [u,v,w]  =  r^pP[f (0 ,0) ,g(9 ,0) , 
h(9,0],  Eqs.  (2.1a-c)  can  be  expanded  further  in  the  series 

03 

[u,v,w]  rXpP[f(0,0)  ,g(e,0),h(9,0)] 

p=0,i,l,. . . 


where  f(9,0),  g(9,0),  h(0,0)  and  X  can  be  obtained  for  each  value  of 
p,  for  a  properly  chosen  function  p  (such  as  p  =  sin  0).  Thus,  form¬ 
ing  a  more  general  represention  of  the  displacement  field  which  is 
easily  accessible  to  numerical  methods  for  its  solution  as  shown  in 
Section  2.2  B-D.  Indeed,  the  above  equation  is  not  the  most  general 
expression  because  there  exists  the  possibility  of  some  other  re¬ 
presentation  for  the  displacement  which  might  be  intuitively  obtained. 
However  complicated,  general,  or  exact  series  representation  one  may 
choose,  the  most  interesting  and  practical  term  in  the  proposed  series 
is  that  whose  stress  field  dominates  over  all  other  possible  terms 


present  in  the  complete  solution.  Since  displacements  behave  like 
and  stresses  like  r^  ^  as  r-»0,  the  smallest  value  for  X  must  be 
sought. 

2.3  Methods  for  the  Eigenvalue  Search 

The  problem  is  to  find  the  smallest  eigenvalue  X  such  that 
Re(X)  >  -Jf.  Again,  X  is  limited  to  Re(X)  >  for  the  strain  energy 
to  be  bounded  near  the  point  of  singularity.  Various  methods  might 
be  available  for  the  eigenvalue  search  connected  with  Eq.  (2.14) 

[3,28].  Two  methods  are  presented,  but  only  the  most  efficient  one  is 
chosen  to  solve  the  problem. 

A.)  Quadratic  polynomial  on  the  eigenvalue  problem. 

It  is  interesting  to  note  that  the  stiffness  coefficients  of 
k„(x),  Eq.  (2.13)  are  quadratic  polynomials  in  X,  see  Eq.  (2.11). 
Hence,  the  matrix  k„(X)  may  be  written  as 

[klj]=k  =  a+bX  +  £X2  (2.22) 

where  £,  _b,  and  £  are  real  square  matrices  independent  of  X  and  of  size 
(M  x  M),  M  being  the  number  of  nodes.  So,  Eq.  (2.14),  k  X  =  0,  takes 
the  form: 


a  X  +  X  b  X 


+  X2  c  X 


=  0 


(2.23) 


where  X  is  the  column  matrix  of  X^.  This  expression  is  a  nonsymmetxic 
quadratic  eigenvalue  problem  on  X 


In  an  iterative  eigenvalue  search, k^  would  have  to  be  repeatedly 
evaluated  for  various  X-values.  Obviously,  Eq.  (2.22)  allows  reduction 
in  the  number  of  computations  needed  to  obtain  k^,  since  it  sufficies 
to  determine  three  matrices  of  sizes  (M  x  M)  independent  of  X,  and 
then  evaluate  k  given  X  till  Eq.  (2.1k)  is  satisfied. 

Furthermore,  it  is  useful  to  observe  that  complex  eigenvalues  of 
Eq.  (2.1k)  or  (2.22)  can  occur  only  in  conjugate  pairs,  because 


k(x)x  =  gxl  X  *  k(X)  X 


(2.2k) 


so  that  if  k(X)X  *  0,  then  also 


k(X)X  *  0 


(2.25) 


where  a  superposed  bar  denotes  a  complex  conjugate. 

B. )  Conversion  to  non -homogeneous  system  of  equations. 

The  method  for  the  eigenvalue  search  described  in  this  section 
has  been  used  with  success  in  connection  with  other  problems  which  lead 
to  the  equations  of  the  type  of  Eq.  (2.1k),  for  the  real  case  in 
Ref.  [J,k]  and  for  the  complex  case  [6].  The  method  used  herein  will 
be  explained  in  complex  arithmetic  which  can  easily  be  converted  to 
real  arithmetic  by  ignoring  all  the  imaginary  components. 

Eq.  (2.1k)  represents  a  large  system  of  M  homogeneous  linear 
algebraic  equations  for  the  values  X^  which  belong  to  the  nodes 
j  *  1,2,... ,M. 


k6 


M 

Y  kij(X^)Xj  1  3,  i  =  1,2,...  ,M.  (2.26) 

i-1 


For  a  given  value  of  Poisson's  ratio  v,  the  root  X  and  the  correspond¬ 
ing  eigenvector  X^(n  =  1,2,...,M),  is  evaluated  by  the  following 
technique : 

First,  the  matrix  k„(X,v)  is  calculated  for  one  chosen  value  of 

X  =  XR  +  iXj,  where  XR  =  Re(X)  and  X^  =  Im(X).  Then,  the  equation 

belonging  to  one  of  the  unknowns  X  ,  e.g.,  for  the  term  X  ,  at  the 

n  m 

surface  node,  is  deleted  from  the  matrix  k^,  and  stored  separately. 

This  equation  is  then  replaced  by  the  equation  Xm  =  (1,1)  which  makes 

the  equation  system  non -homogeneous. 

This  matrix,  for  the  new  system  of  equations,  is  non-singular, 

because  X  is  a  simple  root  when  Im(X)  f  0.  Thus,  all  Xn,(n  *  1,2,...,M) 

can  be  solved  by  converting  standard  library  subroutines  for  banded 

real  matrices  to  complex  arithmetic,  see  Appendix  C.  In  the  case  of 

real  roots,  the  new  matrix  is  normally  non-singular. 

Once  the  unknowns,  X  ,  are  solved,  the  right-hand  side,  0  .  of 

n  to’ 

the  original  m-th  equation  is  evaluated.  The  quantity  0  may  be  re- 

garded  as  a  function  of  X„  and  XT,  i.e. : 

R  L 


M 

y  k  x 

L.  mj  j 
j 


(2.27) 


After  a  second  quantity  of  is  evaluated  for  another  value  of  X, 
chosen  at  the  beginning  of  this  procedure,  the  iterative  "regula 
falsi"  method  is  applied  to  \  and  \  separately,  in  order  to  find 
the  value  of  X  such  that  will  be  zero,  i.e.,  Re(Qffi)  =  IniCQ^)  *  0, 
thus,  satisfying  Eq.  (2.14).  Computationally,  Eq.  (2.14)  will 
obviously  not  be  satisfied  exactly.  In  the  program  the  iterative 
"regula  falsi"  method  is  used  until  the  difference  of  two  consecutive 
values  of  X  is  of  order  0(10  ^) ,  and  to  whose  last  value  yielded 
|QJ  -  0(lCf6)*. 

Accuracy  and  convergence  can  be  improved  when  the  m-th  equation 

is  chosen  such  that  )  |k  .j  is  the  largest  of  all  )  |k..|,i  =  1,...,K, 

j  J 

[6].  The  search  for  this  value  was  not  necessary  here  since  all  of 

them  were  found  to  be  of  the  same  order.  Nonetheless,  the  convergence 

of  the  iteration  method  is  sometimes  quite  slow.  To  obtain  a  good 

initial  guess,  X  must  be  scanned  in  small  steps,  usually  of  about  0.05 

and  0.005  for  Che  real  and  imaginary  parts  of  X,  respectively. 

Fig.  (5-6)  gives  an  indication  of  the  sharply  varying  slope  of  Q  v.s. 

m 

X  for  an  example  whose  solution  is  real.  Thus,  much  care  is  required 
to  avoid  missing  the  smallest  root  and  to  keep  the  computational  time 
to  a  minimum. 

The  root  search  subroutine  can  be  generalized  further  for  various 
cases  which  will  be  used  in  later  chapters.  Note  that  the  search 
for  root  X  may  be  geometrically  interpreted  as  the  intersection  of  the 
line  of  solution  for  a  constant  number  of  finite  elements,  N,  with 
the  vertical  line  v  •  constant,  see  Eq.  (2.15)  and  Fig.  (3-5).  For 

*  -2 
The  difference  in  initial  guesses  for  X  is  in  the  order  0(h  )  and  their 

corresponding  Q  values  usually  range  in  the  order  0(Vr)  to  O(h^). 


regions  when  this  curve,  X  v.s.  v,  turns  sharply  upwards  (or  downwards) 
the  subroutine  converges  poorly  or  not  at  all,  i.e.,  either  the  inter¬ 
sections  occur  at  very  small  angles  or  no  intersection  seems  to  exist. 
To  circumvent  this  difficulty,  Eq.  (2.  16)  may  be  considered  as  an 
eigenvalue  problem  for  v  at  a  fixed  X-  Then,  the  solution  represents 
an  intersection  of  the  line  of  contant  N  with  the  horizontal  line 
X  *  constant. 

Similarly,  when  the  problem  is  to  find  orientation  angles  8  and  y 
for  their  crack  front  edges  and  planes,  or  notches  with  opening  a  and 
orientation  8  and  y,  for  which  the  values  of  X  and  v  are  desired,  as 


fixed  then  k^  will  be  a  function  of  these  angles: 


kij  *  kij(X,v;  tt»8,v) 


(2.28) 


and  Eq.  (2.27)  may  be  considered  as  an  eigenvalue  problem  for  a,B, 
and  y. 

In  the  case  of  complex  roots  the  search  method  is  more  complicated 


than  that  for  the  real  roots.  First,  one  must  scan  a  region  of  com- 

.  2  .  .  2  -1 


plex  X  *  (XR,Xj).  At  each  point  the  quantity  [ReCQ^)  +  Im(Qm>  ] 
is  computed  and  a  plot  is  constructed.  Once  a  peak  is  noticed  in  the 
plot,  then  the  "regula  falsi"  method  is  performed  in  the  following 
manner : 


Two  values  for  \  ,  chosen  inside  the  region  of  the  peak,  are 
fixed.  For  each  of  the  two  X„,  XT  is  iterated  with  respect 
to  Im^)  and  when  it  converges  Re(Qm)  is  stored.  The  third 
XR  value  is  iterated  with  respect  to  the  two  previous  Re(Qm). 


k9 

Again,  it  is  fixed  and  Xj  is  iterated  untill  X^  converges. 

And  so  on.  Once  X,,  has  converged,  a  final  iteration  is 
R 

performed  on  X^.  Thus,  a  complex  root  is  approximately 

computed. 

C. )  Eigenvalue  convergence. 

Interesting  results  can  be  obtained  when  the  convergence  rate  of 
the  numerical  eigenvalues,  computed  for  different  number  of  grid  sub¬ 
divisions,  is  studied  carefully.  Indeed,  convergence  studies  can  at 
times  be  questionable,  especially  when  sophisticated  manipulations 
are  needed  to  obtain  numerical  values.  However,  its  use  is  justified 
when  a  convergence  method  constantly  agrees  with  exact,  or  nearly 
exact,  solutions  to  the  same  problem  based  on  a  completely  independent 
approach . 

For  example,  it  is  well  known  that  the  ordinary  finite  element 

method  exhibits  quadratic  convergence  (or  it  has  an  error  of  the  order 
2 

0(h  ))j  provided  that  there  are  no  singularities  within  the  domain 
[26].  That  is,  functions  f,  g,  and  h  and  their  gradients  are  non¬ 
singular,  p  ■  £,  see  Section  (II-2).  So,  the  convergence  of  the 
eigenvalue  X  should  also  be  quadratic.  When  singularities  are  pre¬ 
sent,  e.g. ,  p  *  0,  then  the  convergence  is  less  than  quadratic  and 
the  error  E  is  of  the  order  0(hm) ,  where  h  is  the  size  of  the  element 
and  m  the  convergence  rate.  Then,  noting  that  (^denoting  proportionality) 

h2  ~  l/N  (2.29) 


and 
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E  ~  k/N  (2.30) 

where  N  is  the  number  of  finite  elements  and  k  a  constant,  possibly 
dependent  on  0 

Q=2v/(1-2v)  (2.31) 

These  relations  should  hold  accurately  when  N  is  sufficiently  large. 
Hence , 


E=XN-Xe~k/N  (2.32) 

where,  X^  is  the  computed  value  using  N  finite  elements  and  Xg  the 
exact  solution.  Then,  Eq.  (2-30)  can  be  written  as 

log  E  3  log  k  +  m  log  ./ft  (2.35) 

If  a  quadratic  convergence  is  present,  expression  (2.3I)  must 
exhibit  a  straight  line  of  slope  m  3  2  for  sufficiently  large  N  when 
log  E  is  plotted  versus  log  /ft.  Otherwise,  |m|  <  2.  This  observation 
can  be  used  to  advantage  in  extrapolating  the  convergence  pattern 
and  estimating  the  results  for  N  -»  ®,  or  h  -»  0. 

D. )  An  extrapolation  technique. 

Many  extrapolation  techniques  exist  in  the  literature,  most 
referring  to  particular  problems  [23,2k]  based  on  the  hm  extrapolation. 
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This  idea  was  first  suggested  by  Richardson  [29]  and  a  fuller  treat¬ 
ment  was  later  given  by  Richardson  and  Grant  [8].  The  latter  deviced 
an  extrapolation  formula,  better  known  as  the  "deferred  approach  to 
the  limit",  in  which  h  represents  the  average  size  of  the  interval 
divisions.  This  method  may  be  extended  here  to  suit  the  problems 
in  question. 

If  is  the  solution  at  the  end  of  an  interval  obtained  by  using 
h  »  h^  m  l//Yl  and  \ ^  is  the  solution  at  the  end  of  the  same  interval 
using  the  same  formula  but  with  h  ■  hg  ~  l//l^,  the  extrapolation 


X 


extrap. 


(2-54) 


gives  an  improved  approximation  over  the  linear  extrapolation  (m  *  1) 
provided  that: 

i. )  the  total  round-off  error  is  negligible 
ii. )  both  h  are  small  enough  for  the  error  to  be  proportional 
to  hm,  i.e.  E  *  0(hm). 

When  N  is  too  large,  there  is  a  danger  that  round-off  error  will  build¬ 
up  to  substantial  proportions.  Thus  far,  this  error  has  not  yet  been 
detected,  even  with  largest  system  of  equations  used  here:  975  simultane¬ 
ous  equations  belonging  to  288  elements.  Hence,  condition  (i.)  may  be 
considered  to  hold  for  the  large  computers  available  today. 


52 


the  convergence  race  m,  In  condition  (ii.)  is  a  seldomly  known 
value,  yet,  it  is  of  most  importance  for  extrapolation  studies.  Then, 
the  question  arises  for  its  value  which  would  apply  to  these  problems. 
Its  assessment  will  new  be  analyzed. 

Extending  the  Richardson's  hm  extrapolation  formula,  Eq.  (2.32), 

which  is  based  on  the  two  conditions  mentioned  earlier,  the  convergence 

pattern  can  be  exploited  further  to  greatly  improve  the  accuracy  of 

the  results  with  the  additional  provision  that  the  grids  for  various 

subdivisions  are  all  similar  and  generated  according  to  the  same  rule 

imposed  at  the  beginning  of  the  problem.  Let  X  be  the  root  obtained 

N 

when  N  number  of  finite  elements  is  used.  The  following  extrapolation 
technique  is  proposed. 

n/2 

Plot  the  values  of  root  XN  versus  1000/n  '  for  various  chosen 
values  of  n.  The  convergence  rate  m  *■  n  which  gives  the  best  straight- 
line  fit,  as  indicated  by  the  least  sum  of  absolute  deviations,  is 
selected.  Then,  a  regression  line  is  passed  to  obtain  the  extra¬ 
polated  value  as  N  -»  ®,  or  h  -»  0. 

Obviously,  this  technique  must  work  if  the  assumption  that  the 
error  is  of  the  order  0(hm)  *  0(N  m^)  holds.  Note  that  all  eigenvalues 
XN  are  included,  not  just  two,  as  in  Eq.  (2.32).  Hence,  the  value  m, 
obtained  in  such  manner,  is  the  effective  convergence  rate. 

This  technique  can  be  interpreted  in  the  following  manner: 

Let  *  1,2, ...,K,  be  the  roots  obtained  by  using  the  finite 

element  method,  e.g. ,  N  “  128,72,32,18;  then  i  *  1,2,3,^.  Let 
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-  h"  -  1000/Nn//2 


(2-55)  ' 


be  Che  X-coordinate  of  Che  X^-root  for  a  convergence  race  n,  where  n 
can  vary  continuously.  Construct  a  rectangular  coordinate  system 
X-x,  as  shown  in  Fig.  3-4.  The  data  points  in  Fig.  3 .k  are  from 
a  typical  example  Co  be  studied  in  the  next  chapter. 

The  best-fit  straight-line  X  *  a  +  bx  through  the  data  points 
is  determined  by  the  minimum  value  of  the  sum  of  squares  of  the 
diviations 

K  K 

S(n)  [Xi  -  a(n)  -  b(n)Xi(n)]  (2.36) 

i*l  i»l 

where  a  =  a(n)  and  b  «•  b(n)  are  the  coefficients  of  the  straight-line. 
Henceforth,  the  subscript  i  will  be  dropped  because  the  summations  are 
understood  to  be  over  all  K  data  points.  For  S  to  be  a  minimum  for 
a  particular  n-value: 


as 

3a 


(2.37) 


which  yields 


b 


'Zxilx 


(2.38a) 


5^ 


i(Xxt-bTxt> 


(2.38b) 


The  optimum  value  of  n,  for  which  S(n)  is  a  total  minimum  is  obtained 
by 


~  S(n)  *  0,  at  n  “  m 


(2-39) 


from  which  m,  the  effective  convergence  rate,  is  intuitively  chosen. 
Clearly,  the  extrapolated  value  will  be  a(m).  The  proofs  that: 


E  **  0(hm)  ■  S(n)  *  0,  at  n  *  m 
on 


(2.40) 


or  vice  versa 

~  S(n)  »  0,  at  n  »  a  ■  E  *  0(hm)  (2.4i) 

on 

are  beyond  the  scope  of  this  work,  if  such  theorem  does  indeed  exist. 

Unfortunately,  the  literature  on  this  technique  is  not  available. 
Its  servicability  can  only  be  supported  if  it  concurs  with  known 
solutions,  as  it  will  be  shown  to  be  true  in  the  succeeding  chapter. 

E. )  Consent. 

The  application  or  a  regular  finite  element  method  to  the  vari¬ 
ational  equation,  Eq.  (1.14),  is  straight-forward,  with  the  only  im¬ 
plication  being  computational  errors.  The  technique  used  to  search 


for  the  smallest  complex  eigenvalue  has  been  proved  to  work  success¬ 
fully  in  [6].  An  extrapolation  technique  based  on  convergence 
patterns  is  proposed.  A  rigorous  proof  that  such  technique  must  work 
is  not  given. 


CHAPTER  III 


NUMERICAL  SOLUTIONS 


3- 1  Introduction 

The  finite  element  method  developed  on  the  (8-0) -plane  and 
applied  to  the  variational  equation  derived  on  the  unit  sphere  is  of 
general  applicability.  This  method  can  now  be  used  to  obtain  the 
solution  for  a  crack  whose  .front  edge  and  plane  inclination  angles  are 
of  arbitrary  values,  see  Fig-  3.18,  page  103;  a  notch  of  arbitrary 
opening  and  orientation,  see  Fig.  3.26,  page  101,  and  for  the  solution 
of  a  crack  or  notch  with  two  dissimilar  materials. 

The  analytical  solution  for  a  Mode  I  crack  whose  front  edge  and 
plane  are  normal  to  the  surface,  see  Fig.  3.1,  page  68,  has  been  ob¬ 
tained  by  Benthem  [9]  and  Kawai,  Fugitaui,  and  Kumagai  [10].  Signi¬ 
ficant  advances,  which  led  to  highly  accurate  analytical  solutions, 
have  recently  been  made  in  potential  theory  problems  by  Morrison  and 
Lewis  [11]  and  by  Keer  and  Parihar  [12].  The  former  authors  succeeded 
in  obtaining  a  tractable  differential'  equation  by  virtue  of  using 
special  coordinates  (conical  coordinates)  suited  for  the  particular 
problem  of  charge  singularities.  Keer  and  Parihar's  method,  utilizing 
spherical  coordinates,  appears  to  have  broader  application  and  involves 
the  use  of  Green's  functions  to  formulate  the  problem  in  terms  of  a 
singular  integral.  The  crucial  step  is  to  differentiate  this  integral 
equation  to  get  rid  of  a  constant  right-hand  side  and  obtain  an  eigen¬ 
value  problem,  which  is  then  solved  numerically  by  Erdogan  and  Gupta's 
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method,  and  thus  obtaining  the  solution  for  crack  corner  in  an 
infinite  elastic  space  in  Mode  I  opening;  also  obtained  numerically  by 
Bazant  [3],  who  used  finite  difference  methods  as  an  approach  to  the 
problem.  Parihar  and  Keer  have  extended  their  very  effective, 
original  and  elegant  method  to  the  same  problem  for  Modes  II  and  III 
singularities  which  is  irreducible  to  potential  theory  [14a] .  They  have 
also  obtained  the  solution  for  shear  on  a  rigid  corner  stamp  on  a  semi¬ 
infinite  elastic  body  for  which  the  solution  is  complex  [14b].  These 
solutions  and  those  of  plane  problems  provide  valuable  check  cases 
for  the  accuracy  and  correctness  of  the  present  method.  In  a  more 
recent  private  communication*  Benthem  has  obtained  numerical  solutions 
not  yet  published  of  an  arbitrary  crack  using  finite  difference  methods 
applied  to  the  differential  equations  of  equilibrium.  His  solutions 
agree  reasonably  well  with  the  results  to  be  presented. 

In  the  progress  of  this  work  certain  limitations  to  the  finite 
element  method  have  been  found.  The  obvious  one  is  that  for  which  the 
Poisson  ratio  is  close  to  0.5  and  the  term  Q  =  2v/l-2v  increases  with¬ 
out  bounds,  and  for  which  it  was  noted  that  Modes  II  and  III  are  more 
susceptible  than  Mode  I.  Also,  when  the  angle  of  inclination  8  for  the 
crack  front  edge  is  close  to  0  or  rr,  see  Figs.  3. 10  and  3. 11,  page  84, 
numerical  inaccuracies  were  seen;  because  when  these  domains  are  mapped 
in  the  (9-0) -plane  they  are  distorted  considerably  and  one  would  need 
to  increase  the  number  of  finite  elements  until  the  domains  are  rea¬ 
sonably  represented.  A  final  limitation  is  that  whenever  the  eigen¬ 
value  is  real  and  larger  than  unity  in  the  interation  routine,  \  will 
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converge  Co  exactly  unity,  the  reason  being  that  rotational  effects 
will  dominate,  i.e.,  X  *  1. 

3-2  Check  Cases 

As  a  first  step,  the  program  is  checked  for  its  correctness  and 
accuracy.  Various  simple  cases  of  known  solution,  usually  given  in 
terms  of  displacements  in  the  Cartesian  coordinate  system,  are  trans¬ 
formed  to  the  spherical  coordinate  system  [7,  page  37 ]•  This  is  done 
by  letting  the  y-axis  coincide  with  the  crack  plane  0  ■  0,  tt  at  9  =*  tt/2; 
the  z-axis  coincide  with  the  crack  front  edge,  0  *  0;  and  the  x-axis 
being  perpendicular  to  both,  y-  and  z-axis,  i.e.,  at  0  *  +  tt/2, 

0  «•  it/2,  see  Fig.  (3*1)  page  68.  Then,  the  following  transformation 
is  allowed: 


- 

1 

i 

1 

sin  0 

sin  0 

sin  0 

COS  0 

cos  0 

| 

1 

v  r 

\ 

= 

cos  0 

sin  0 

cos  0 

COS  0 

-sin  0 

u  f 

y 

i 

w 

L 

cos  0 

-sin  0 

o 

J 

1  \ 

where  (ux,uy,u2)  is  the  Cartesian  displacement  field  of  the  known 
solution. 

The  spherical  displacement  field,  (u,v,w),  in  the  domain  0  £  0  * 
rr/2,  0  S  0  S  rr  is  checked  for  i)  continuity,  ii)  existence  of  at  most 
first  order  derivatives,  and  iii)  boundary  condition  requirements. 
Then,  the  field  is  substituted  into  the  program  by  calculating 
the  displacements  at  each  nodal  point,  i.e.,  obtaining  X^,  Eq.  (2.14). 


A 


The  stiffness  matrix  is  subsequently  computed;  and  Eq.  (2.lk)  must 
be  approximately  satisfied.  As  an  error  indicator,  the  right-hand 
sides  for  all  i  were  compared  to  the  sum  of  their  absolute  terms  as 
indicated  by  the  condition 

M  M 

I  I  kijXj|/ I  lkijXjl  <  10'4  (3'2) 

i=l  i=l 

at  all  nodal  points  i,  i  =  1,2, for  a  mesh  of  only  32  elements. 
The  elementary  solutions  for  the  various  special  cases  considered  here 
were  first  analyzed  for  their  dependence  on  X  and  Eqs.  (2.1a-c), 
in  order  to  obtain  the  X  and  p  values. 

A. )  Rigid  body  rotations. 

The  three  body  rotations  allowed  imply  that  X  *  1  and  p  =  0.  For 
example,  the  rotation  about  the  z-axis,  8=0,  implies  that 

u=v=0,  w  ■  r  sin  9  y  (3- 3) 

for  which  X  *  1,  p  *  0,  H(0,0)  =  sin  9;  and  F(0,0)  *  G(6,0)  *  0.  Note 
that  it  is  also  possible  to  have  X  *  1,  p  =  1,  and  h(8,0)  =  1  for  this 
particular  example.  Table  3*1  shows  the  print-out  of  Eqs.  (2. lV)  and 
(3.2)  using  Eq.  (3-3)- 

B.  )  Homogeneous  strain  field. 

The  only  homogeneous  stress  field  that  will  satisfy  the  free 
surface  conditions  is  that  which  belongs  to  Oyy  *  1  (or  constant),  i.e. 


Table  3.1:  Nuraberical  values  at  nodal  points  for  rigid  body  rotation 


h(S,$)  *  -sin  $  cos  $ 


(3.6c) 


RHS  FORCES  IN  REDIRECTION 
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Table  3.2:  Numerical  values  at  nodal  points  to  homogeneous  strain  field 
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would  apply.  Table  3-2  shows  the  print-out  of  Eqs.  (2.110  and  (3*2) 
using  Eqs.  (3-6a-c). 

C. )  Plane-strain  solutions. 

The  solutions  for  near-tip  plane-strain  fields  nay  be  found  in 
references  [  ]  for  opening  I  mode,  shear  mode  II,  and  antiplane 

mode  III.  Obviously  the  antiplane  mode  field  cannot  satisfy  all  stress 
boundary  conditions  at  the  surface  6  =  tt/2,  or  the  nodes  which  belong 
to  the  body  surface  in  Eq.  (2.110.  In  this  case  only  the  fulfillment 
of  the  equilibrium  equations  for  the  Interior  nodes  was  checked.  As 
an  example,  the  displacement  field  for  mode  I  opening  is  [12]: 


Qfr  [sin  0(A  +  B)] 

(3-7a) 

C/r  [cos  0 (A  +  B)] 

(3-7b) 

C/r  [A-B] 

(3-7c) 

where  C  -  -Kj/2  (HvVe1;  El  =  e/(1-v2);  v1  =  v/(l-v) 

1  2 

A  *  [2(l-v  )  -  cos  a]  sin  a  sin  0 

(3-8) 

1  2 

B  *  [1  -  2v  +  sin  a]  cos  a  sin  0 


a  *  (0  -  tt)/2 
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Table  3.3:  Numerical  values  at  nodal  points  for  near-tip  plane  strain  field 
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Then,  X  ■  i  and  p  *  0.  These  values  also  hold  for  the  two  remaining 
modes  and  need  not  be  discussed  further.  A  print-out  of  Eqs.  (2.1k) 
and  (3*2)  is  shown  in  Table  3-3  using  Eqs.  (3.7a-b)  and  (3*8). 

D.  )  Checks  on  smaller  domains. 

For  the  sake  of  accuracy  for  the  method,  cases  B  and  C  were 
rerun  for  domains  of  small  notches,  but  still  containing  32  elements. 
Table  3.4  shows  the  print-out  for  the  example  given  in  Section  C  for 
a  notch  of  boundaries  0  £  9  £  tr/l6,  15rr/l6  <i  0  SS  tt-  As  expected, 

Che  accuracy  increases  inside  the  domain  but  not  at  the  boundaries, 
since  the  boundaries  of  the  actual  problem  are  not  those  of  a  notch. 

E.  )  Comment. 

Note  that,  if  Eq.  (2.l4)  is  satisfied  computationally,  i.e., 
its  rlghthand  sides  are  small,  then  alternatively,  the  variational 
equation,  Eq>  (1.13),  must  be  satisfied  exactly.  In  all  check  cases 
studied  above  substitution  of  Eqs.  (3.7a-c),  (3-5a-c),  and  (3.3)  into 
Eq.  (1.13)  yielded  zero  after  long  hand  algebraic  manipulations 

3*3  Crack  Plane  and  Front  Edge  Normal  to  Surface 

The  finite  element  computer  program  derived  in  Chapter  II  and 
outlined  in  Appendix  B  is  now  applied.  The  first  problem  is  that  of 
a  crack  whose  plane  and  front  edge  are  normal  to  the  halfspace,  as 
depicted  in  Fig.  3.1.  Recently,  Benthem  [9]  and  Kawai,  Fujitani,  and 
Kumagai  [10,30]  presented  analytical  solutions  for  this  problem  but 
only  for  Mode  I  opening.  A  comparison  of  their  results  with  the  ones 
obtained  in  this  chapter  is  made. 
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A- )  Symmetric  opening.  Mode  I;  p  * 

To  analyze  the  field  near  the  terminal  point  0,  Fig.  3*1,  in 
Mode  I  opening  for  the  problem  just  presented  it  is  sufficient  to 
consider  only  half  the  domain,  because  there  exists  symmetry  with 
respect  to  0  =  tt-  The  new  domain  will  continue  to  be  rectangular  in 
the  (0-0)-plane  with  boundaries  0  £  0  £  tt/2,  0  £  0  £  tt,  Fig.  3*2,  or, 
as  indicated  by  the  domain  enclosed  by  the  slashed  lines  in  Fig.  3*1. 

The  stress  boundary  conditions  on  the  crack  surface  (0  =  0)  and  on  the 
half -space  surface  (0  =tt/2)  are  automatically  satisfied  by  the  finite 
element  method.  The  boundary  conditions  of  0*0  (the  pole,  top  side 
of  the  (0-0)-domain  in  Fig.  3*2),  are  irrelevant  and  none  have  been 
imposed. 

The  boundary  conditions  on  the  symmetry  plane  (0*tt)  must  properly 

reflect  the  symmetries  of  displacements  and  stresses  with  respect  to 

0* it.  Therefore,  for  the  symmetric  crack  (Mode  I)  opening,  one  must 

impose  for  all  nodes  at  0=rr  the  condition  w  =  0,  i.e.,  h  =  0,  Eqs.  (1.6c) 

and  (2.1c).  The  symmetry  conditions  for  stresses,  namely  o.  =  (J,Q  =  0, 

0  r  0U 

will  be  also  be  automatically  satisfied  by  the  finite  element  method  as 
natural  boundary  conditions.  Thus,  these  considerations  ensure  a 
statically  determinate  support  for  the  body  and  at  the  same  time  properly 
reflect  the  symmetry  properties. 

From  the  work  previously  done  on  potential-related  problems  [3],  it 
was  expected  that  the  displacement  field  should  exhibit  a  behavior  of 
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Fig.  3- 


Fig.  3 


Orthogonal  crack.  Spherical  coordinate 
system  at  termination  of  crack  front  edge 
00’ at  body  surface,  pointO.  (The  unit 
sphere  is  shown  only  to  visualize  the  co¬ 
ordinate;  the  body  is  semi-inf inite) . 


.2:  Finite  element  grids  used  for  orthogonal 
crack.  Domain  O’ACO'  from  Fig.  3-1 
visualized  in  the  (0-®) -plane. 


v/n 

18 

32 

72 

128 

Benthem  [9] 

0.0 

0.565263 

0.537891 

0.517198 

0.50973 

0.50 

0.15 

0. 648844 

0.611861 

0.582320 

0.570591 

0.5164 

0.30 

0.756772 

0. 704-681 

0.662787 

0.645832 

0.5477 

0.40 

no  conver. 

0.82639e 

0.756209 

0.721745 

0.5868 

x/n 

18 

52 

72 

128 

0.905 

0.390605 

0.423676 

0.453383 

0.466796 

Table  3*5 :  Numerical  results.  Eigenvalues  for 

orthogonal  crack  using  N  finite  elements; 
Mode  I,  p  = 
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u(r,6,0)  -  r^sinP9  y  (6,0)  (3*8b) 

In  Ref.  [3],  p  was  suitably  taken  as  p  ■  i.  This  choice  was  motivated 
by  the  fact  that  the  term  pp,  Eq.  (3.86)  is  dominant  at  finite  r  and  as 
0-»O.  Unfortunately,  the  literature  related  to  the  method  of  solution 
used  here  to  treat  elasticity  problems  is  non-existent.  Therefore, 
expecting  similar  behaviors,  p  was  chosen  to  be  p  B  J  as  a  first  at¬ 
tempt  to  solve  the  problem. 

Symmetric  opening  is  acquired  in  the  program  by  forcing  any  of  the 
nodes  belonging  to  the  crack  surface  (0  *  0)  in  the  0-direction. 

Table  3*5  gives  the  numerical  results  of  \  for  various  values  of 
Poisson's  ratio  v.  For  values  of  v  which  exceed  0.4,  the  root  search 
surbroutine  converged  very  poorly,  or  not  at  all.  For  these  cases,  X. 
was  fixed  and  v  was  considered  the  root,  as  explained  in  Section  2.3  B. 
Some  results  are  given  at  the  bottom  of  Table  3-5 •  However,  when  v 
becomes  very  close  to  0-5 ,  the  present  formulation  breaks  down,  because 
the  value  Q  *  2/(1  -  2v) ,  Eq.  (1.8),  increases  without  bounds.  A 
special  program  would  have  to  be  written  for  v  close  to  0.5  and  for 
incompressible  materials,  v  ■  0.5- 

Note  that  for  the  case  of  Poison's  ratio  v  ■  0  the  computed  value 
of  the  root  for  the  finest  grid  used  (128  elements,  439  simultaneous 
equations),  was  0.50973-  The  exact  solution  is  known  to  be  0.5  [9]- 
Thus,  the  computed  value  is  still  within  1-9#  error.  Closer  estimates 
for  the  exact  solution  with  these  values,  Table  3-5,  can  be  gotten 
with  the  extrapolation  technique  explained  in  Section  2-C.  Eq.  (2.33 ) 
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should  hold  for  m  ■  -2,  since  p  =  §  will  not  introduce  gradient  singu¬ 
larity  near  8  -*  0,  [26].  The  plot  of  log  E  versus  log /N,  Eq.  (2*33)» 
is  shown  in  Fig.  3. 3  for  the  case  v“0,  where  the  exact  solution  is 
known,  Xg  *  0.5*  Indeed,  the  plot  is  a  straight  line  with  a  slope 
indicating  in  =  -2.0.  Thus,  for  v  =  0,  the  present  formulation,  p  ■ 
seems  to  follcw  a  systematic  pattern  of  quadratic  convergence. 

This  observation  can  be  used  to  advantage  in  extrapolating  the 
convergence  pattern  and  estimating  the  results  for  N  ■*  ®,  h  -»  0;  even 
for  v  >  0,  where  no  error  analysis  can  be  made,  since  no  exact  solution 
is  available.  Thus,  a  plot  of  X  versus  1000/N,  i.e.,  m  =  2,  is 
constructed  in  Fig.  3*4.  Again,  for  quadratic  convergence  these  plots 
should  be  straight  lines  for  sufficiently  large  N.  According  to 
Fig.  3*4  this  seems  indeed  to  be  true.  Therefore,  regression  lines 
(straight  lines)  are  extended  to  obtain  estimates  of  the  values  as 
N  -*  ®,  i.e.,  estimates  of  the  exact  solution,  as  shown  in  Fig.  3>4. 

The  extrapolated  values,  along  with  the  numerical  results  of  Table  3.4, 
are  shown  In  Fig.  3-5  and  ate  compared  with  Benthem's  solution  [9]. 

Note,  however,  that  for  the  case  v  “  0,  the  extrapolation  point,  N  -»  ®, 
falls  on  0.5  +  0.002,  X  ■  0.5  being  the  exact  solution. 

The  fact  that  estimates,  N  -»  ®,  significantly  deviate  from 
Benthem’s  results  [9],  as  shown  in  Fig-  3-5,  can  be  attributed  to  the 
case  p  i. )  The  solution  presents  eigenvalues  which  are  in  the  order 

of  those  obtained  by  Benthem,  but  unfortunately  for  the  case  p  =  the 
exact  or  numerical  solution  is  unlikely  to  be  available  for  comparison 
purposes.  ii. )  In  light  of  the  results  to  be  presented  in  the  sub¬ 
sequent  section,  this  solution  is  correct  within  1  %.  iii. )  From  (ii.), 


N  =  18 


Fig-  3’J:  Determination  of  the  rate  of  convergence 
with  increasing  number  of  elements.  Use 
of  Eq.  (2.33);  Mode  1,  p  ■ 


J 


Fig-  3-4:  Extrapolation  of  numerical  results  to 
infinite  number  of  elements,  using 
Eq.  (2.31);  Mode  I,  p  *  § . 


(3x6) 
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Fig*  3*5:  Singularity  exponent  X  for  various 
values  of  Poisson  ratio*.  Mode  I. 
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this  solution  assures  the  existence  of  an  infinite  enumerable  eigen¬ 
values  for  the  problem,  e.g.,  other  solutions  can  be  gotten  for 
p  *  0,^,1,...,.  Hence,  Fig.  3-5 >  Is  not  a  complete  solution,  unless 
V  m  0. 


B.)  Symmetric  opening,  Mode  I,  p  *  0. 

The  author  is  obliged  to  John  P.  Benthem,  Professor  at  Delft 
University  of  Technology,  for  clarifying,  in  a  private  communication, 
the  value  for  p  from  the  implications  of  his  analytical  solution. 

The  choice  p  *  \  in  previous  computations.  Fig.  3*5 >  was  in¬ 
appropriate  for  the  complete  solution,  because  a  restriction  is  pre¬ 
scribed  to  the  displacements,  similar  to  that  of  a  generalized  Fourier 
series  which  would  represent  the  displacements,  thus,  limiting  their 
complete  and  natural  dependence  on  the  angle  0  for  which  the  smallest 
eigenvalue  should  exist.  Let  (rp)P,  Eqs.  (2.1a-c)  be  the  tern  with 
the  lowest  exponent  in  the  field  near  the  singularity  line,  (crack  front 
edge  p  “8  *  0).  Indeed,  p  *  If  is  the  lowest  p  corresponding  to  the 
deformed  states  for  crack  front  singularity,  but  where  the  displacement 
field  behaves  like: 


u,v,w  m  9^  ;  9  -*  0  ,  0  <  r  <  »  . 


and  stresses  (displacement  gradients),  like: 
o  ~  9~^  ;  9  -»  0  ,  0<r<®. 


(3-9) 


(3-10) 
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However,  in  the  neighborhood  of  the  singularity  line  (crack  front 
edge,  6>-  0)  one  may  have  values  of  p  *  0,  *4,  1,  ...  etc.,  as  was 
proven  in  Section  2.2k,  Eq.  (2.6).  In  otherwords,  Eq.  (3.9)  does 
not  mean  that  there  are  no  displacement  fields  starting  with  the 
stronger  exponent  p  *  0, 

u,v,w  ~  8°  ;  0  -»  0,  0  <  r  <  ®  .  (5*11) 

The  exponent  p  ■  0  does  not  give  rise  to  stresses 

a~61;8-*0,0<r<®.  (3-12) 

In  otherwords,  the  term  (p)P  does  not  cause  any  singulatiry  as  p  -»  0, 
or  0  -*  0,  at  a  finite  fixed  r.  However,  this  may  cause  functions 
F (0,0),  and  H(0,0)  to  have  gradient  singularity  of  the  type  0p+<*,  or 
0^,  as  r  -*  0,  where  q  >  -1.  This  singularity  would  be  more  severe 
than  the  singularity  0”^  associated  with  the  planar  near  tip  field. 

q 

That  terms  of  9  ,  q  >  -1  as  r  -*  o  should  indeed  be  present  is  indicated 
by  Benthem's  solution  [9 ]•  This  will  still  satisfy  the  restriction 
that  along  the  crack  front  edge  the  strain  energy  must  remain  finite, 
i.e.  ,  the  behavior  0**  of  the  stresses  along  the  crack  front  edge  be 
such  that  q  >  -1. 

Therefore,  all  finite  element  solutions  were  rerun  with  the 
exponent  p  *  0.  The  numerical  results  are  given  in  Table  3.6,  and  are 
compared  with  Benthem's  results.  Again,  as  expected,  the  error  de- 
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creases  considerably  by  using  larger  number  of  elements.  Also,  the 
more  Poisson's  ratio  v  deviates  from  zero  the  larger  the  error,  due 
to  the  value  of  Q  ■  2/(1  -  2y),  Eq.  (1.17a).  Obviously,  Q  will  Induce 
numerical  inaccuracies  which  are  well  known  to  occur  in  all  numerical 
methods  of  this  type. 

It  was  fortunate  that  for  the  results  obtained  in  the  previous 
section,  p  ■  J,  the  convergence  rate  m  was  known  "a  priori'  and 
proved  grafically  to  be  quadratic  by  employing  Eq.  (2.31)*  For  the 
case  p  “  0,  for  which 

f(0,0)  =>  F(0,0),  g(0,0)  -  G(0,0),  h(O,0)  =  H(0,0)  (3-13) 

m 

while  p  a  J  is  also  present,  the  convergence  rate  must  be  less  than 
quadratic.  But,  since  the  exact,  or  nearly  exact  solution  is  available 
[9],  Eq.  (2.31)  can  again  be  used  grafically  to  find  the  value  of  m, 
as  shown  in  Fig.  3.6.  The  extrapolation  points  (or  regression  lines 
N  -*  ®) ,  are  shown  in  Fig.  3*7.  These  points  are  then  compared  with 
Benthem's  solution  in  Fig.  3.8,  showing  both  solution  coinciding  with 
each  other  within  a  0.002  deviation. 

Because  the  gradient  of  F(0,$),  G(0,$)  and  H(0,$)  might  tend  to 
infinity  as  0  -*■  0,  it  seems  appropriate  to  refine  the  grid  step  A6  as 
0  decreases.  Irregular  rectangular  net  works  in  which  A<|>  was  con¬ 
stant  and  in  which  A0  was  refined  so  as  to  keep  A0  roughly  equal 
(sin  0)A$,  have  been  tried,  using  same  numbers  of  subdivisions  in 
both  0  and  $  directions,  as  shown  in  Fig.  3.9.  Although  the 
numerical  results  for  the  same  maximum  size  element  (regular  grids) 


Table  3.6:  Numerical  Results.  Eigenvalues  for  orthogonal  crack 
using  N  finite  elements;  Mode  I,  p  -  0. 


1000  /  N 

Extrapolation  of  numerical  results  frost  Fig.  }.6 
to  K  *  «.  Mode  1.  p  »  0. 


Fig.  3.8:  Numerical  results  for  orthogonal  crack. 
Mode  I,  p  *  0. 
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were  better  than  those  for  the  refined  (irregular  grids), 
the  extrapolated  values  for  the  refined  grids  were  no  better  than 
those  for  uniform  subdivisions.  So,  non-uniform  subdivisions  of  the 
meridians  would  be  Ineffective. 

Recently,  Kawai,  Fujitani,  and  Kumagai  [10]  also  presented 
analytical  solutions  for  the  same  problem  (orthogonal  crack).  They 
obtain  three  roots  for  all  Poison  ratios  which  disagree  with  Benthem 
[9]  as  well  as  the  present  work.  For  example,  the  smallest  root  [10] 
for  v  ■  0* 3  is  approximately  X  “  0*3 •  The  insert  of  Fig.  3.6  shows 
the  value  of  Q,  which  must  be  zero  in  the  search  of  the  eigenvalue, 

Eq.  (2.26).  The  curve  that  Q  versus  X  traces  is  smooth  and  continuous 
Therefore,  no  eigenvalue  near  X  “  0. 3  for  Q  -  0  could  have  been 
missed. 

Furthermore,  the  program  was  checked  against  the  analytical  and 
numerical  solutions  for  a  sharp  corner  of  angle  2a  on  the  crack  front 
edge  of  a  planar  crack  whose  complement  is  the  wedge-shaped  punch 
of  angle  2q  within  an  infinite  elastic  solid,  see  Fig.  3* 10.  The 
solution  for  this  synmetric  opening  (Mode  I)  of  such  a  crack  was 
given  in  [3],  where  a  finite  difference  solution  was  based 
on  a  reduction  to  potential  theory.  Very  accurate  solutions,  by 
means  of  singular  integral  equations,  have  recently  been  obtained 
by  Keer  and  Parihar  [13].  Both  solutions  [3,13]  have  found  that  for 
symmetric  opening,  the  eigenvalue  is  independent  of  Poison  ratio 
for  a  fixed  angle  2a. 


and  Extrapolation  to 


Fig.  3.9:  Extrapolation  of  numerical  results  for  orthogonal  crack 
with  grid  refinement  for  Poisson  ratio  v  =  0. 15. 


3* 10:  Crack  corner  in  an  infinite  elastic 
space.  (The  unit  sphere  is  not  the 
body  surface*  it  is  used  to  visualize 
the  spherical  coordinates j  the  body 
is  infinite.  ) 


3* 11:  Finite  element  grid  used  for  the  crack 
corner  in  an  infinite  elastic  space. 
Domain  O'ACO'  from  Fig.  3-  10  visualized 
in  the  (8 -0) -plane. 
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Because  of  symmetry  About  $  ■  jt  and  0  *  8,  one  needs  only  to 
consider  one  quarter  of  the  unit  sphere  as  shown  by  the  dashed  lines 
of  Fig.  3* 10.  The  bottom  side  of  the  domain,  0^  *  a/2,  is  a  great 
circle  around  the  unit  sphere,  which  in  the  (9-0)-plane  is  given  by 
the  equation 


6  •  *  arctan  (tan  g/cos  #);  if  0  <  0,  then  0  «-  8  +  tt  (3*1^) 

where  g  ■  a/2.  The  0-coordinates  of  the  nodal  points  on  the  r-th  curved 
row.  Fig.  3* 11]  are  calculated  as 

0  -  0b(r-l)/(rb-l)  (3.15) 

where  rfc  is  the  number  of  the  last  row,  0  *  0^,  and  r  ■  1  corresponds 
to  0  ■  0.  Eq.  (3-15)  describes  a  uniform  subdivision  of  each  meridian, 
as  shown  in  Fig.  3-  11* 

The  stress  and  displacement  conditions  at  0  ■  tt  have  already 
been  discussed  in  the  previous  section.  The  displacement  conditions 
at  the  bottom  boundary  (0^),  must  be  replaced  by  displacement  boundary 
conditions  of  symmetry: 

v  cos  15  -  w  sin  n  *  0,  at  0  *  0^  (3*  16) 

where  n  is  the  angle  that  the  bottom  boundary  makes  with  the  d-axis 

in  the  (@-0)-plane,  i. e.  ,  the  normal  displacement  of  the  bottom  nodal 

points  is  zero.  The  stress  boundary  conditions  a  *  or  -  0  at 

00  0r 


8' 


9  *  9^  will  be  automatically  satisfied  by  the  finite  element 
method. 

Table  3*  7  gives  the  numerical  results  for  the  cases  of  v  “  0 
and  0.3  for  a  ■  n/^»  chosen  for  examples.  Using  Eq.  (2.31)  the 
values  of  m  are  graphically  calculated  in  Fig.  3>  12.  Finally,  the 
extrapolation  point  are  obtained  in  Fig.  3* 13  >  and  compared  to  the 
values  X  ■  0.296,  calculated  in  [3]  and  \  ■  0. 2966  in  [13].  Indeed, 
the  eigenvalues  are  independent  of  Poisson  ratio  for  a  fixed  a.  The 
fact  that  all  three  values  (X  s  0. 296  obtained  by  different  researchers 
using  independent  methods  of  solution),  are  the  same,  further  confirms 
that  the  present  solution  is  correct. 

Now,  the  question  arises  for  the  value  of  the  convergence  rate 
m  when  no  exact  solution  is  available  and  hence  Eq.  (2. 31)  cannot  be 
used  as  before.  For  such  cases,  the  extension  of  Richardson's  hm 
technique,  described  in  Section  2. 3D,  was  found  to  work  exceptionally 
well.  All  numerical  results  were  run  in  a  simple  subroutine,  which 
is  included  in  the  program,  Appendix  B.  The  m-values  are  given  in 
the  preceding  Tables,  next  to  their  corresponding  numerical  values. 

The  extrapolation  values  for  all  cases  thus  far  studied  came  within 
0.  Ulf>  error,  as  shown  in  the  tables.  Again,  reaffirming  the  present 
method  of  solution  and  justifying  the  use  of  the  extrapolation  tech¬ 
nique  proposed  in  Section  2. 3D.  For  the  sake  of  brevity,  whenever 
an  extrapolated  value  is  mentioned  herein  it  will  refer  to  this 
technique,  unless  otherwise  specified. 


log  ( X  - 


-  log  ( X  - 
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Note  that  the  extrapolation  value  in  the  subroutine  corresponds 
to  the  first  coefficient  of  the  straight-line  with  optimum  slope  m, 
i.  e.  ,  the  value  of  the  regression  line  at  h  *  0,  N  ». 

A  final  check  was  performed  on  the  orthogonal  crack  in  Fig.  J>.  1, 
for  Mode  I  opening.  It  consisted  of  including  the  entire  (0-0) -domain 
of  the  unit,  sphere,  i.  e. ,  0  s  0  £  rr/2,  0  s  0  s  2rr;  where  no  symmetry 
considerations  need  be  made.  One  node  may  be  fixed  in  the  0-dlrectlon 
to  prevent  rotation  and  thus  implement  a  statically  determinate  sup¬ 
port  for  the  body.  However,  since  rotation  implies  X  *  1,  see  Section 
5-2A,  the  support  is  normally  not  necessary  unless  X  “  !•  This  also 
means  that  in  this  finite  element  method  it  is  sufficient  to  impose 
only  one  force  at  a  nodal  point  to  achieve  the  mode  required,  and 
not  two  forces  of  opposite  direction  applied  at  two  opposite  nodal 
points. 

The  numerical  results  are  given  in  Table  5>  8  along  with  the 
convergence  rate  m  and  the  extrapolated  value  for  the  case  v  *  0.  15, 
but  excluding  the  value  for  N  *  18.  If  this  last  value  were  to  be 
included  one  would  obtain  m  *  2. 6  and  X  *  0.5267.  However, 

note  that  N  *  18  in  the  domain  0  s  9  £  ri/2,  0  £  0  £  2-^  is  a  very 
coarse  mesh,  thus,  Inducing  an  error  which  would  not  be  of  the  order 
0(hm),  m  <  2.  Since  the  convergence  rate  is  to  be  limited  by 
m  £  2,  as  mentioned  earlier  in  Section  2,3C  the  value  for  N  *  18 
has  to  be  excluded,  even  though  the  extrapolated  value  using  all  four 
points  is  within  a  2$  error.  And  as  a  rule  of  hand,  so  will  future 
values  for  N  ■  18  when  the  entire  domain  is  Included.  For  such  cases, 
the  number  of  finite  elements  will  be  raised  to  200  and  288,  where 


i 


one  would  again  expect  a  convergence  rate  less  than  quadratic  and 
more  accurate  results. 

C.  )  Antisymmetric  openings,  Mode  II  and  III. 

For  antisymmetric  crack  openings  the  question  of  proper  anti¬ 
symmetric  conditions  at  0  ■  tt  and  at  the  free  surface  9  *  ti /2  are 
more  complicated  than  for  the  symmetric  opening.  It  appears  that 
Modes  II  and  III  cannot  exist  seperately  at  the  surface  point  (which 
was  first  suggested  by  Professor  L.  M.  Keer  of  Northwestern  University 
in  an  uncontested  $5-  00  bet).  Indeed,  it  is  impossible  to  imagine 
conditions  of  zero  stress  state  at  the  half-space  surface  (0  *  tt/2)  bo 
be  satisfied  by  a  displacement  field  which  would  exhibit  either 
Mode  II  or  Mode  III  antlsynsnetry.  The  finite  element  claculations 
confirmed  this  also-  i.  e.  ,  when  the  full  domain  0c(O,ti/2),  0e(O,2fr) 
was  used  and  Mode  II  antisymmetric  displacements  were  forced  in  two 
symmetrically  opposite  nodes  at  the  crack  surface  (n  ■  1  at  0  ■  tt/2 
and  n  *  -1  at  0  ■  5tt/2,  both  at  0  ■  n/2),  the  v  displacements  at 
0  *  tt/2  were  found  to  be  nonzero  and  exhibit  perfect  antisymmetry 
about  0  ■  tt;  which  is  characteristic  of  Mode  III.  Furthermore,  a 
surface  nodal  displacement  was  forced  such  that  Mode  III  opening 
would  be  obtained,  i.  e. ,  the  v  displacement  at  the  crack  surface. 
However,  the  eigenvalue  slowly  converged  to  the  same  eigenvalue  when 
Mode  II  opening  was  forced.  Thus,  the  antisymmetric  Modes  II  and  III 
are  always  combined  at  the  surface  point. 

Therefore,  one  may  impose  at  0  ■  tt»  the  symmetry  plane,  either 
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Mode  II-type  condition:  u  sin  0  -  v  cos  0*0,  or  (3*17) 

Mode  Ill-type  condition:,  u  co  0  -  v  sin  6*0,  (3*18) 

or  any  linear  combination  of  these  two  conditions,  among  which  the 
simplest  choice  is 

u(9,tt)  “  v(0,tt)  “  0.  (3-19) 

The  antisymmetry  condition  for  stress  in  both  modes  is  a  *  0  at 

00 

00  *  tt,  which  is  again  automatically  satisfied  by  the  finite  element 
method  as  a  natural  boundary  condition. 

The  singularity  exponent  \  in  either  case  is  the  same,  and 
because  it  belongs  to  a  combination  of  two  modes,  X  Is  a  double  root. 

Table  3-9  gives  the  numerical  results  when  Eq.  (3-17)  is  used 
and  Table  3* 10  gives  the  numerical  results  when  the  full  domain 
8e(0»Tl/2),  0c(O»2rr)  is  used  where  no  symmetry  analysis  is  made.  The 
extrapolated  values  in  Table  3* 10  differ  significantly  from  those 
of  Table  3- 9»  Hence,  at  this  point  it  was  necessary  to  increase  the 
number  of  finite  elements  to  288  for  the  domain  0e(O,2|r),  as  shown  in 
Table  3* 10  for  the  case  v  ■  0. 3-  This  case  yielded  the  same  extra¬ 
polated  value  as  that  obtained  from  the  domain  0e(O,n)*  Fig.  3. 14 
shows  the  numerical  results  of  Table  3* 10  with  the  extrapolated 
values  of  Table  3-  10.  Note  that  for  the  case  v  *  0  the  eigenvalue 


32 


72 


128 


m 


Benthem 


v/n  18 


X 


extrap. 


0.15  0.789175  0.639575  0.572169  0.549754  1-99  0.5234  0.5164 


Table  3*  8:  Numerical  results.  Eigenvalues  for 

orthogonal  crack  with  full  body*  0  £  g  £  n/2, 
0  £  0  £ 2tt;  using  N  finite  elements. 

Mode  I,  p  «  0. 


v/N 

18 

32 

72 

128 

m 

^•extrap. 

0.0 

0.612712 

0. 564782 

0.529639 

0.517067 

1.926 

0.50001 

0. 15 

0.551476 

0-500448 

0. 463557 

0.450516 

1. 966 

0.43533 

0.3 

0.521529 

0. 466533 

0.426113 

0. 411730 

1.922 

0.  40207 

0.4 

0.530050 

0. 465491 

0.415043 

0.398513 

1.860 

0.39591 

Table  3*9:  Numerical  results.  Eigenvalues  for 

orthogonal  crack;  O£0£ti/2,  O£0£tt;  using 
N  finite  elements.  Modes  II  and  III. 


v/N 

32 

72 

128 

200 

288 

m  x 

Aextrap. 

0.0 

0.698366 

0.597145 

0.555897 

1.580  0.48432 

0.  1 

0.653992 

0.545166 

0.504840 

1.842  0.44715 

0. 15 

0.639572 

0. 528586 

0.  487784 

1.862  0.43030 

0.3 

0.616785 

0.  498913 

0. 454053 

0. 432559 

0.420677 

1.880  0.40202 

0.4 

0.635112 

0.506881 

0. 452256 

1.842  0.34790 

Table  3* 10:  Numerical  results.  Eigenvalues  for  orthogonal 

crack  with  full  body;  O£0£jt/2,  O£0s2tt;  using 
N  finite  elements.  Modes  II  and  III. 
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is  exactly  \  *  0. 5,  which  was  again  expected.  No  solution  seems  to 
exist  in  the  literature  concerning  Modes  II  and  III  for  comparison 
in  Fig.  3.  14. 

The  antisymmetric  opening.  Modes  II  and  III,  was  also  checked 
against  the  known  analytical  solution  of  a  crack  corner  in  an 
infinite  elastic  space,  Fig.  3* 10,  page  84  solved  by  Keer  and 
Parihar  [14a],  They  found  that  the  solution  for  these  modes  to  be 
irreducible  to  potential  theory  and  to  depend  on  Poisson  ratio  v 
For  example,  for  a  crack  corner  of  angle  a  *  tt/4  the  eigenvalues 
X  ®  0.2966  and  O.3285  were  obtained  for  Poisson's  ratio  v  *  0.  0 
and  0.25,  respectively.  Table  3. 10a  gives  the  numerical  results 
for  both  cases.  Figs.  3* 15  shows  the  convergence  pattern  using 
Eq.  2.33  and  Fig.  3-16  the  extrapolated  values  using  the  convergence 
rate  m  obtained  graphically  from  Fig.  3*  15-  The  values  for  m  and 
^extrap  08 1*1®  ^.39  are  also  given  in  Table  3*  10a. 

3.  4  Crack  Propagating  at  the  Surface 

From  the  practical  point  of  view,  the  case  of  a  propagating 
crack  is  of  main  interest.  There  exist  certain  physical  restrictions 
for  the  solution  of  a  propagating  crack  which  can  be  derived  from 
energy  considerations.  For  cracks  that  do  not  propagate,  the  only 
restrictions  are  that  the  strain  energy  within  a  small  sphere  about 
point  0,  as  well  as  the  strain  energy  per  unit  length  of  edge  within 
a  small  cylinder  whose  axis  coincides  with  the  crack  front  edge  00', 
Fig.  3* 1»  he  integrable.  Let  the  strain  energy  by  denoted  by  Eq, 


then 


finite  elements.  Modes  II  and  III, 
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Eo  ■  J  %  ®1J  «lj  dv- 

v 


(3.20) 


The  above  expression  can  be  simplified  by  asymptotic  analysis. 
Noting  that  near  the  surface  terminal  point  0,  the  displacements 
behave  like: 


u 


i 


~  r 


X 


(3.21a) 


where  t*  denotes  proportionality.  Hence, 

au/ar  ~  r*'1 

(3.21b) 

«ij  ~  rx 

(3.21c) 

_  rX-l 

aiJ  ~  r 

(3. 2 Id) 

dv  M  r^dr 

(3. 21e) 

yields. 

I  -  r^+1 
o 

(3.22) 

For  the  strain  energy  to  be  integrable,  or  bounded,  as  r  0,  re¬ 
quires 
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Re(x)  >  -%  ,  stationary  crack.  (3-23) 

As  the  crack  propagates,  energy  flows  into  all  points  of  the 
crack  front  edge  and  is  consumed  by  the  process  of  separation,  1.  e.  , 
creation  of  crack  surfaces.  The  energy  flux  near  the  points  of  the 
crack  front  edge  may  generally  have  two  components:  (a)  The  flux 
E^  which  is  parallel  to  the  edge  and  flows  into  any  point  on  the 
crack  front  edge,  including  the  surface  point  0.  must  be  zero 

because  the  trace  of  the  surface  point  0  as  it  moves  is  a  line, 
and  a  line  can  be  associated  only  with  a  negligible  amount  of  ad¬ 
ditional  surface  energy.  (b)  The  flux  Eg  of  energy  into  the  moving 
crack  front  edge  per  unit  length  of  edge  must  be  finite  and  non¬ 
zero  because  the  surface  energy  y  is  finite  and  non-zero. 

The  first  condition  (a)  requires  that 

Ei  *  II  aij(3uj/3*)  dn  =  o  (3-24) 

where  is  the  cartesian  stress  tensor;  are  the  cartesian  dis¬ 
placements*  x  is  the  coordinate  in  the  direction  of  the  crack  ex¬ 
tension;  and  q  is  a  surface  of  a  sufficiently  small  sphere  centered 

2 

at  point  0.  Noting  that  df}  «  r  sin  0  d0d<$,  and  from  Eqs.  (3.21),  it 

2\ 

follows  that  Ej  w  r  ,  and  for  E^  to  be  zero  as  r  -»  0,  it  is  necessary 
that 


Re(x)  >  0 


» 


(3-25) 


101 


which  is  a  weak  condition  on  a  propagating  crack. 
The  second  condition  (b)  requires  that 


\  -  J*  fo^/ax)  r^0 


(3-26) 


where  (r^,0)  is  a  polar  coordinate  system  in  a  plane  normal  to  the 
crack  £ront  edge*  L  is  a  circle  of  radius  r^  in  this  plane  centered 
around  the  edge;  x  is  the  direction  of  crack  propagation;  are 
the  cartesian  displacements;  and  <j^  is  the  cartesian  stress  tensor. 
Also,  the  energy  flux  Eg  may  in  general  be  expressed  by  Rice's 
J-integral  [23,31]  for  linear  elastic  behavior: 


T  9Ui 

h 58  JL  &  au  *ijdy  *  °ij  \  sr  rid«) 


(3-27) 


in  which  dy  -  r^  sin  0  d0  and  a/dx  -  cos  0(d/9r^)  -  (sin  <t>/r^) 

(b/S0).  On  physical  grounds,  the  flux  Eg  must  obviously  be  positive, 
non-zero  and  finite  at  all  points  near  the  surface  terminal  point  0 
of  the  crack  front  edge.  Furthermore,  the  flux  Eg  may  be  expected 
to  be  constant  along  the  crack  edge,  assuming  that  the  energy  needed 
for  the  creation  of  new  surface  is  the  same  along  these  points. 
However,  this  last  requirement  may  be  simplified  by  the  asymptotic 
deductions.  When  Eqs.  (3*21)  are  substituted  into  Eq.  (3.27),  it 
follows  that  (**  denotes  proportionality). 


(3-28) 


Thus,  £or  Eg  to  be  bounded  and  non-zero  as  -»  0,  it  is  necessary 
that  Re(2\  -  1)  »  0,  or 

Re(X)  *  propagating  crack  (2.3 9) 

This  condition  oust  be  satisfied  for  the  terminal  surface  point  of  a 
crack  that  propagates,  but  not  for  a  stationary  crack,  as  it  is  well 
known  [23,31]. 

According  to  Eq.  (3.29),  a  crack  which  propagates,  or  for 
which  propagation  is  imminent,  must  exhibit  X  “  £  (the  exponent  being 
assumed  to  be  real  if  there  are  no  two  dissimilar  materials).  By 
far,  this  case  is  of  the  greatest  practical  interest.  Therefore, 
a  meaningful  question  is  to  ask  whether  there  exist  inclinations  g 
of  the  crack  front  edge  Fig.  3* 17  and  y  of  the  crack  plane  Fig.  3*  18  , 
for  which  the  eigenvalue  X  *  %  is  attained.  For  the  orthogonal 
crack  edge  (g  *  n/2,  y  *  0;  Fig.  3-l)»  propagation  is  obviously  pos¬ 
sible  only  if  v  *  0. 

Tables  3.  11  and  3. 12  give  the  numerical  results  for  the  symmetric 
(Mode  I)  and  antisymmetric  (Modes  II  and  III)  cracks  whose  plane  is 
normal  to  the  surface  (y  ■  0)  and  whose  edge  Inclination  angle  0 
varies  for  various  values.  The  extrapolated  results  are  plotted  in 
Fig.  3-20. 

It  is  interesting  to  note  that  the  solutions  presented  in  Fig.  3*20 
agree  with  the  common  sense  that  as  the  crack  "size"  defined  by  the 
edge  inclination  angle  g  decreases,  1.  e. ,  there  is  more  material 
that  is  not  cracked,  the  eigenvalue  X  increases;  and  as  the  crack 
"size"  Increases,  1.  e. ,  there  is  more  material  that  is  cracked,  the 
eigenvalue  X  decreases.  In  other  words,  the  stress  singularity 
exponent  X-l  increases  (weaker  singularity)  as  g  decreases  and 


1T;  Spherical  coordinate  system  at  termination 
of  crack  front  edge  00*  at  body  surface. 
Inclined  edge.  (The  unit  sphere  is  shown 
only  to  visualize  the  coordinates;  the  body 
is  semi-infinite). 


. 18:  Spherical  coordinate  system  at  termination 

of  crack  front  edge  00’  at  body  surface.  In¬ 
clined  edge  and  inclined  crack  plane.  (The 
unit  sphere  is  shown  only  to  visualize  the  co 
ordinates-  the  body  is  semi-infinite). 
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inclination  angle  p. 


Table  5. 11:  Numerical  results.  Eigenvalues  for  crack  whose  plane  is 
normal  to  the  surface  and  whose  edge  inclination  angle 


Table  5- 13:  Numerical  results.  Eigenvalues  for  crack  front  edge  angle 
g  of  a  propagating  crack,  v  ■  0. 3-  Mode  I. 


decreases  (stronger  singularity)  as  £  increases.  Also  note  that 
the  lines  in  Fig.  J.20  are  not  straight  lines  and  they  seem  to  ap¬ 
proach  X  ■  1  for  g  -♦  0,  i.  e.  ,  there  is  no  crack  and  one  would  ex¬ 
pect  rotational  effects  to  take  place,  see  Section  $.2  Aj  and  they 
seem  to  approach  \  *  0  for  g  -*  tf,  i.  e.  ,  the  semi-infinite  body  is 
completely  cut  in  half  and  one  would  expect  rigid  body  translations 
to  take  place. 

In  order  to  substantiate  the  accuracy  of  these  results,  the 
approach  to  the  eigenvalue  problem  was  modified  by  treating  the 
stiffness  matrix  as  a  function  of  angle  $  rather  than  X>  i. e. ,  X  was 
fixed  to  0.5  and  Eq.  (2.27)  was  treated  as: 

M 

^kij(g)Xj=0,  J«1,2,...,M  (3-29) 

j-1 

The  eigenvalue  search  routine  based  on  the  Newton  method  was  easily 
converted  to  search  for  0  instead  of  X>  This  alternate  method  was 
tried  for  the  case  v  =  0*  3  In  Mode  I  opening  and  drawn  separately  in 
Fig.  3.21(a).  The  numerical  results  obtained  by  this  method  and  based 
on  up  to  288  elements  are  given  in  Table  3.13.  The  convergence 

pattern  of  g  versus  N  is  shown  in  Fig.  3- 21(b)  along  with  the  extra¬ 
polated  value  (N  -*  ®)  which  yields  the  same  value  as  that  obtained 
grafically  in  Fig.  3.21(a).  Therefore,  the  values  of  g  for  each 
value  of  Poisson's  ratio  for  which  X  =  0. 5  from  Fig.  3*20  are  drawn 


0.60 


V 

Fig.  J. 22:  Dependence  of  crack  front  edge  angles  p 
of  a  propagating  crack  upon  Poisson's 
ratio  v.  (a)  Mode  I*  (b)  Modes  II  and 
III.  (Normal  crack  plane). 
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The  physical  meaning  of  the  edge  inclination  angle  g  for  which 
X  *  0.5,  is  that  for  a  propagating  crack  the  symmetric  opening. 

Mode  I,  gives  an  obtuse  angle  (p  >  tt/2),  i.  e.  ,  the  surface  pont  0 
trails  behind  the  interior  crack  edgej  and  the  antisymmetric  opening, 
Modes  II  and  III,  gives  an  acute  angle  (g  <  n/2),  i.  e. ,  the  surface 
point  0  moves  ahead  of  the  interior  crack  edge.  The  fact  that  they 
are  different  has  an  important  physical  consequence:  At  the  terminal 
point  of  a  crack  whose  plane  is  normal  to  the  surface  a  combined 
mode  propagation  is  impossible,  i.  e.  ,  the  crack  would  assume  such  a 
shape  that  its  surface  terminal  point  propagates  either  with  a  sym¬ 
metric  opening.  Mode  I,  or  with  an  anitsymmetric  opening,  Mode  II  and 
III,  but  not  both  combined. 

In  view  of  this  result,  it  is  natural  to  ask  whether  there  exist 
an  inclination  angle  y  of  the  crack  plane  for  which  the  p-values  for 
the  symmetric  and  antisymmetric  exitation  of  a  propagating  crack 
(X  s  0-5)  would  coincide,  see  Fig.  3.  18.  However,  the  numerical  results 
given  in  Tables  3-1^  and  3-15  for  the  case  v  *  0. 3  and  drawn  in 
Fig.  3’ 23  indicate  that  this  never  occurs,  and  as  the  crack  plane  becomes 
inclined  (d  ±0),  the  3  -  values  for  Re(A)  m  V2  vary  as  a  function  of  6. 

In  these  cases  it  is  no  longer  possible  to  distinguish  between  symmetric 
(Mode  I)  and  antisymmetric  (Modes  II  and  III)  openings,  for  there  is  no 
geometrical  symmetry.  For  each  of  the  two  S-values,  there  exists  at 
point  0  a  certain  limiting  ratio  t^:  Kji  K3  of  the  stress  intensity 
factors  for  Modes  I,  II  and  III  and  no  other  ratios  are  possible.  So, 
for  cracks  of  inclined  plane,  the  propagation  of  the  surface  point  takes 
place  always  in  a  combination  of  all  three  modes.  Conversely,  for  a 
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given  ratio  1^:  K^,  one  can  generally  find  the  angles  B  and  6 
which  must  get  established  at  the  surface  point. 

In  solving  this  problem  one  must  take  into  account  the  entire 
domain  0  s  0  s  g,  y  ^  0  ^  2^  +  y,  because  the  symmetry  is  destroyed 
when  y  >  0.  Fig.  3* 19  shows  a  typical  finite  element  grid  in  the 
(0-0) -plane  corresponding  to  the  domain  enclosed  by  the  dashed  curve 
of  Fig.  3. 18. 

3.5  Experimental  Fracture  Specimens 

Some  recently  obtained  experimental  results  allow  a  check  on  the 
present  numerical  results.  These  are  the  fatigue  loading  fracture 
tests  made  by  P.  D.  Bell  and  W.  J.  Feeney  [15],  to  whom  the  author 
is* obliged  for  making  their  results  available,  and  are  reproduced  in 
Figs.  3-  2h  and  3*  25.  These  photographs  show  the  crack  arrest  marks 
observed  in  fatigue  Mode  I  fracture  tests  of  alluminum  alloy  and 
titanium  alloy  specimens.  The  Poisson  ratios  of  these  materials  are 
(according  to  material  handbooks)  about  0.35  and  0.32,  respectively, 
and  for  which  the  present  solution,  Fig.  3*20  gives  8  «  102°  for  both 
materials.  These  angles  are  plotted  and  compared  in  Figs.  3*2^  and 
3.25.  Comparatively,  the  observed  trend  agrees  with  fhe  numerical 
resutls  in  that  the  surface  point  trails  behind  the  interior  crack 
edge  (i. e. ,  0  >  90°)  rather  than  moving  ahead.  The  numerical  value 
does  not  agree  too  closely  with  the  observed  average,  but  considering 
that  some  small  scale  yielding  and  Inelastic  strain  reversals  occur 
in  the  actual  tests,  and  that  the  plastic  "shear  lip"  phenomenon  can 
along  cause  0  >  90°,  the  comparison  cannot  be  qualified  as  poor. 

One  must  also  realize  the  inevitable  statistical  scatter  of  the 
experiment. 


Fig.  3.24:  View  of  fracture  face  showing  fatigue  arrest 
marks  on  Aluminum  Alloy  2219-T5  8l  (y  -  0. 33 
Mode  I  crack,  sheet  width  6.35  ran  (l/4  in.  ), 
magnification  7* 4  times,  crack  propagates  to 
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The  crack  arrest  marks  indicate  the  line  of  constant  X  *  §  for 
which  propagation  is  established.  The  curves  that  these  marks  trace 
may  be  explained  in  the  following  manner.  Prior  to  propagation,  the 
two-dimensional  theory  of  fracture  mechanics  gives  X  ■  \  for  the 
interior  of  the  specimen  (plane  strain  in  the  middle  regions  and 
plane  stress  near  the  surface)  while  \  >  \  will  hold  at  the  surface 
point.  When  propagation  is  established  the  interior  edge  (x  * 
will  move  perpendicularly  to  the  crack  Independent  of  Poisson  ratio, 
while  the  surface  point  (weaker  singularity,  X  >  £)  will  have  to 
reach  X  *  £  by  allowing  the  crack  edge  to  shift  angles  and  thus  moving 
behind  the  interior,  as  predicted  in  Fig.  3-20. 

Even  though  experimental  results  are  nonexistent  for  Modes  II 
and  III,  a  similar  reasoning  may  be  made.  Again,  the  two-dimensional 
theory  of  fracture  mechanics  gives  X  *  i  for  the  interior  points 
while  X  <  ^  will  hold  at  the  surface  points.  When  propagation  is 
established,  provided  the  propagation  plane  remains  in  the  rame  plane 
as  that  of  the  crack,  the  surface  point  (stronger  singularity  X  <  ^) 
will  have  to  reach  X  *  \  by  allowing  the  crack  edge  to  shift  angles 
and  thus  move  ahead  of  the  interior,  as  predicted  in  Fig.  3.20. 

3-  6  The  Two-Material  Interface 

In  plane  elasticity,  the  singularity  exponent  of  an  interface 
crack  between  two  dissimilar  materials  is  complex.  Consequently, 


the  displacements  in  a  close  enough  neighborhood  of  the  crack  tip 
oscillate  along  the  radial  ray.  This  implies  an  overlap  of  crack 
faces  which  is,  of  course,  physically  impossible  and  is  prevented  by 
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contact  of  crack  surfaces.  Nevertheless,  it  is  generally  believed 
that  the  field  for  complex  \  ir  at  least  applicable  in  not  too  close 
neighborhood  of  the  crack  tip,  well  beyond  the  region  of  oscillations. 
That  this  is  indeed  the  case  for  Mode  I  cracks  has  been  demonstrated, 
by  Comninou  [22].  It  must  be  noted,  though,  that  recently  more 
physically  meaningful  solutions  which  take  into  account  the  contact 
stresses  on  crack  surfaces  have  been  developed  [22],  but  their  adaptation 
is  beyond  the  scope  of  this  program.  Thus,  while  extension  of  these 
developments  to  three  dimensional  singularities  should  be  of  high 
priority,  at  present  we  must  be  content  with  the  less  than  perfect 
oscillating  singularity. 

The  foregoing  solution  applies  without  any  change  to  cases 
where  \  is  complex.  Then,  of  course,  and  X^,  Eq.  (2.15),  must 
be  also  considered  complex  and  the  program  must  be  converted  to 
complex  arithmetic,  which  is  easily  achieved  by  proper  type  declara¬ 
tion  of  FORTRAN  variables.  Some  difficulties  were  caused  by  the  need 
of  an  equation  solving  subroutine  for  complex  banded  nonsymmetric 
matrices.  Such  subroutine  has  not  been  available  in  standard  soft¬ 
ware  packages,  and  so  it  had  to  be  developed,  and  it  is  listed  in 
Appendix  C,  page  144. 

A.  )  Check  cases. 

In  the  first  four  sections  of  this  chapter  it  was  noted  that  the 
convergence  of  the  eigenvalue  as  the  number  of  finite  elements  in¬ 
creased  was  systematic  and  an  extrapolation  technique  was  thus  developed. 
In  the  check  cases  where  the  eigenvalues  are  complex,  the  real  part 
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was  found  to  behave  equally  well,  but  the  complex  part  did  not.  In¬ 
stead,  the  complex  part  of  \  gave  nearly  the  exact  value  for  any 
number  of  finite  elements.  These  check  cases  were  the  well  known 
two-dimensional  interface  crack  solutions  and  the  problem  of  a  rigid 
corner  stamp  of  angle  2g  which  was  solved  by  Parihar  and  Keer  [14b]. 
The  two-dimensional  solution  for  the  exponent  is  given  by 

[32,33] 


*■  “  £  ±  iYx 


(3-30) 


where. 


Y.  *  g  log  [  (—  +  — )/  (•—  +  — 
12 


)] 


(3-31) 


and  *  3  “  for  plane  strain,  (since  generalized  plane  stress 
cannot  be  modeled  by  the  program) . 

For  a  study  case,  materials  whose  Young's  moduli  have  the  ratio 
E^/Eg  *  l/hO  and  whose  Poisson  ratios  have  the  same  value  =*  V2  * 

0.3,  were  chosen.  For  these  values  Eq.  (3. 31)  yields  =  0.  O887. 

The  domain  to  be  considered  must  be  0  s  0  <  rr/2,  0  s  0  s  2jt,  where 
the  elements  in  the  region  0  >  tt  have  a  Young  modulus  forty  time  larger 
than  those  in  the  region  0  <  tt,  but  both  regions  have  the  same  Poisson 
ratio  *  Vo  *  0.3;  end  to  simulate  the  two-dimensional  problem  sup¬ 
ports  perpendicular  to  0  *  rr/2  must  be  placed.  The  numerical  results 
were:  \  -  O.65OV7  +  0.  065^1;  0.57329  ±  0.  08528i;  0.5^88  +  0.  08785i 
for  N  ■  32,  72,  and  128  elements,  respectively.  The  extrapolated 
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value  of  the  real  parts  gives  0.49709  which  is  within  0.4%  error  of 
the  exact  value  of  0.5.  Note  that  the  imaginary  part  of  the  eigen¬ 
value  using  128  elements  is  already  within  1.5%  error.  When  v  «  0 
for  both  materials,  the  plane  strain  solution  [33]  also  applies  for 
the  surface  singularity  with  an  orthogonal  crack  edge;  for  a  2:1 
ratio  of  young  moduli  this  gives  X  =  0.5  ±  0.0535i;  whereas,  the 
program  yielded  Im  (X)  -  0.0514  for  N  *  128. 

For  the  second  study  case,  a  rigid  corner  stamp  of  angle  2g  * 
0.2886tt  on  a  semi-infinite  body  of  Poisson  ratio  v  =  0-3,  were  chosen 
from  the  table  in  Ref. [14b]  where  the  analytical  solution  is  given 
by  0. 2U-Ti+  +  0.  0409 i.  The  numerical  results  for  this  problem  were 
X  *  0.37044  +  0.  04393 i‘,  0.3128  +  0.  045 32 i;  0.28804  +  0.  0450Qi  for 
N  -  32,  72,  128  elements,  respectively.  The  extrapolated  value  of 
the  real  parts  is  0.241,  again  within  0.4$;  and  that  the  imaginary 
part  using  128  elements  is  within  10$.  There  may  be  two  possible 
reasons  for  the  imaginary  part  to  be  in  such  relatively  large  error: 
a.  )  The  representation  of  the  exact  domain  with  finite  elements  is 
not  very  accurate  for  such  angle  g,  which  unfortunately  was  the 
largest  angle  that  Parihar  and  Keer  could  consider,  b.  )  The  ana¬ 
lytical  solution  obtained  by  Parihar  and  Keer  involves  an  approximate 
function  substituting  a  Bessel  function,  which  restricts  them  to 
consider  only  small  angles. 

In  either  case,  the  solutions  obtained  with  the  present  pro¬ 
gram  show  that  numerical  results  can  be  obtained  with  reasonable 


accuracy. 
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B.)  Additional  results 

Further  cases  were  run  using  the  present  program  for  which  no 
solutions  have  been  given  before.  These  were  the  results  for  v  >  0. 
i.e.,  the  plane  strain  solution  is  not  applicable. 

a. )  The  singularity  for  an  orthogonal  crack  edge  of  an  inter¬ 
face  of  two  materials  with  a  Young's  moduli  ratio  of  2:1  and  equal 
Poisson  ratio,  the  program  gives  Im(A)  *  ±  0.0399  (N  *  128)  for 

v  ■  0.05  and  Im  (A)  *  ±  0.006  (N  ■  128)  for  v  **  0.3.  It  can  be 
noted  that  Im  (A)  decreases  with  increasing  v.  For  v  =  0.3  and  a 
30:1  ratio  of  E,  the  program  indicated  Im  (A)  to  be  0  or  almost  0. 

b. )  For  these  cases,  the  program  again  showed  Im  (A)  to  be 
close  and  almost  0.  For  an  interior  crack  plane  of  an  orthogonal 
two-material  interface,  the  program  indicated  that  a  crack  with  a 
front  edge  orthogonal  to  the  two-material  interface  has  A  =  0.545, 
0.521,  and  0.499  for  E  -  rations  1:1,  5:0  and  10:0  for  v  =  0.3  (N  -*•  »)  . 

3.7  The  Notch  Surface  Singularity 

Solutions  for  the  surface  singularity  at  notches  have  not  been 
given  before.  The  present  program  can  readily  handle  notches  with 
higher  accuracy  since  the  material  domain  decreases  with  the  size  of 
the  notch  angle.  Fig.  3.26  shows  the  numerical  results  (N  -*■  ®)  for 
notches  terminating  at  the  surface  with  orthogonal  (8*^/2)  and  sym¬ 
metrical  opening  (\>  =  -  a) . 
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CONCLUSIONS 


The  finite  element  method  in  angular  spherical  coordinates  provides 
a  powerful  general  tehcnique  for  determining  three-dimensional  elastic 
stress  singularities.  The  numerical  results  for  cracks  whose  edge  is 
nopmal  to  the  surface  and  for  crack  corners  in  an  infinite  space  are 
in  close  agreement  with  the  analytical  solutions  of  Benthem,  and  Keer 
and  Parihar,  respectively.  The  front  edge  of  a  propagating  crack 
must  terminate  at  the  surface  point  obliquely.  The  values  of  this 
angle  are  different  for  symmetric  (Mode  I)  and  antisymmetric  (Mode  II 
and  III)  crack  opening;  which  indicates  that  a  combined  mode  propagation 
is  impossible  at  the  surface  point  of  a  crack  whose  plane  is  normal  to 
the  surface.  For  Mode  I,  the  surface  point  trails  behind  the  interior 
of  the  crack;  while  for  Modes  II  and  III,  the  surface  point  moves  ahead 
of  the  interior  of  the  crack.  For  cracks  of  inclined  plane,  the 
propagation  at  the  surface  point  takes  place  in  a  combination  of  all 
three  modes. 

The  numerical  results  for  a  rigid  corner  stamp  on  a  semi-infinite 
space  are  also  in  close  agreement  with  the  complex  analytical  solution 
of  Keer.  Some  numerical  results  of  complex  singularities  are  obtained, 
as  well  as  some  cases  of  notches. 
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APPENDIX  A 


DERIVATION  OF  THE 
VARIATIONAL  EQUATION 


The  combination  of  the  equations  of  equilibrium,  Eqs.  (l.Ja-c) 
with  the  boundary  conditions,  Eqs.  (1. lla-c)  is  given  by  the  vari¬ 
ational  statement  Eq.  (1.  12).  With  the  intention  of  using  the 
finite  element  method,  Eq.  (1.12)  has  to  be  reduced  to  an  equation 
which  involves  no  higher  than  first  order  derivatives  and  which 
automatically  includes  the  boundary  conditions.  In  the  following 
derivation  the  asterisk  *  indicates  terms,  or  term,  which  have  been 
added  and  substracted. 


Substituting  Xr,X0,  and  X^  from  Eqs.  (l. Ta-c)  into  the  surface 
integral  of  Eq.  (1.12)  yields 


JJTTCCQCX  -  1)  '2  +  2XHXF  +  2F  +  G0  +  G  cot  9  +  ifrTe  + 


-  C(X  +  1)G0  -  F00  3-  cot  0[(X  +  1)G  -  Fg]  + 

+  T7~~X  (t-t  F^  -  H  -  \H  . )}  sin  9  {>F 
sin  Q  'sin  0  00  0  0  J 


+  [2\G0  +  2\G  cot  0  +  H0  -  (2XG0  +  2XG  cot  0  + 

+  siiTe  sin  eSF  +  t(O+2)(XF0+2F6  +  Gee+Ge  cot  e  - 
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nr?  V  +  Xt(»+1)G  -  re»  sin  650  +  %  ■ 2G  c«ae 


-  ' 20  cot2e)I%i" e6G  +  tnr?  «+2>^»  +  £r4  + 


+  % +  % cot  6  +  nr?  V  -  *<nr?  F«  - E  - *H)  +  <Hee + 


+  H0  cot  6  ‘  7^1  h  +  3f  G«  ‘  V16B  +  ‘-V  6  + 


+  ^TH  +  ^^'Sf  G*  -  *»  -  <-V*  e  ♦ 


I  --g  T1  I  -g-  c  2  cos  9 
sin2  9  *to«  M  sin20 


G0  -  2XH)fsin  0ftH*]]dSd^ 


(A.  1) 


Arranging  some  terms  and  cancelling  others  (a  light  slash  indicates 
terms  which  cancel),  Eq.  (A. 1)  is  reduced  to: 


JJHttQ(X-l)-B]ftF  +  »  +  Ce  +  C  cot  8  +  H  ) 


+  * 2F  +  g9  +  G  cot  9  +  nr?  V -2*Vg  cot  e+nr?  V 


+  (*-Ge  -  Ge  +  f89)  +  cot  eCxc  -  g  +  f0)  + 


+  JiTe  'TITeV^V  V)sln04F 


+  !Q(>-re+2Fe  +  G0e  +G6c°t  9  '  ^  0  + 


+  2(XFe +  (1 + 1)Fe+see  +ceC0t  6  *  ^  G  +  ^  He* 


+  IST  %  •  Vot  6+Jth  V  +  >^  +  1>G  -  U’e 


2G  cot29  -  nib  % +  2G  cot^ 


+  [-Q  cot  9(XF+2F  +  G0  +  H0  +  G  cot  0) 


+  Q  cot  0(XF  +  2F  +  G0  +  0  +  G  cot  0)]*}sin  9  $>G 


+  {ty^’q  (XF  +  2F  +  G  +  G  cot  0  +  -1.-1  H  ) 
sin  0  0  0  @0  0  sin  0  00 y 


+  TTT  Q  (XF  +  (l  +  l)^  +  Gq  +  G  cot  0  +  — r~r  H  ) 

sm  0  0  0  00  0  sin  0  00' 


+  (H0fl  "  Hecot  9  +  - 1—  H  +  -T~r  G  -  -S20-  G  ) 

99  9  sin  0  sin  6  ^  sin29  « 


F^  +  X(X  -  1)H  +  cot  0  (H„  -  H  cot  0  +  —  :  G  ] 
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+  COt  e(H0 


-  H  cot  0  + 


1 

sin  0 


V  +  <2  *75* 

v  sin  0 


2 — V"  +  2*)R  “ 

sin  0 


-  2H*  - 


2  - 
sin  0  G00  +  2XR3sin  0  gH  ] j  d0d0 


(A.  2) 


The  terms  multiplied  by  sin  0  gG  are  arranged  to  form: 


[Q(\Fg  +  2F  +  G„-  +  GgCot  0  -  1 —  G  +  1  -  H  -  ---  0-  H  )  + 

00  00  0  sin2 Q  sin  6  9<$  sin  0  V  + 

+  2(Gee  +  V  +  cot  eCQ(XF+2F+Ge+nr?H0  +  G  cot  0)+ 


+  2tGe  +  F>1  +  iiST  <V  V“  0  +  UT5  s«>  -  2'°'  etrrk h„  + 


'■sin  0  0 


+  G  cot  0  +  F)+\(\+i)g  +  2F  -  Q  cot  0()F  +  2F  +  G  + 

o  0 


+  iiTe  H0  +  G  cot  9)  +  >.F0  -  2( — i—  -  SSfy~  -  l*)G-2G*]si„  esc 

sin  0  sin  0 


(A.  3) 


Placing  Eq.  (A.  3)  back  into  Eq.  (A. 2)  and  rearranging  terms: 


Jit™  '  1)  -2KXF+2F  +  G0+G  cot  9+J^-g  H0)  +2\(\  +  2)F 


+  ^G0  ‘  G0  +  F00  ^  +  cot  0 (XG  -  G  +  F0 ) 
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+  T  F  +  XH  -  H  )lsin  0  6F 

sin  0  sm  Q  00  0  0  J 


+  {-Q  cot  0(XF  +  2F  +  Gg  +  G  cot  0  +  H  )  +  XFg 


-  2  cot  0(,s1^-q  +  F  cot  0  +  F)  +2(XG  -  G  +  FQ)  +  X(0  -  l)G 


+  Q(xre+2F9  +  G89  +  O9cot  e--^-0+niTH6(#.^a-H?>>  + 


+  2«!B9+F0)  +  cot  8[Q(XF+2F+G0+G  cot  H  )  + 


+  2<Vr)1  +  nh'W"  6+1iTe  "00' 'Sitl  0  86 


+  {cot  e(ae-H  cot  e+TiTe  G0)  +2(1^Fe*XH-H) 


”  F  +  X(X  ■  l)H  +  -  H.cot  0  + - H  +-r^  Gfl<4  -  GJ 


sin  0  0 


.96-..e - stn20  "'•‘"a  W  sin20  *' 


+  cot  0(H0-H  cot  6+  G#) 


+  -/Q— ■  (\F  +  2F  ,  +  G.  +  G  cot  9 +  ;■  H  ) 

sin  0  v  0  0  00  0  sin  0  00 


— ~ — -  H  +  G  cot  0  +  F  )]sin  0  JH  T1  dod0 

sin  0  vsin  0  00  0  0 ' i  u  w 


(AA) 


which  can  be  written  as 
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Jit-  e(X-l)-2](XF+2F  +  G6+G  cot  8+ jjj-j  H  ) +2x*ta  e(X+2)F 
+  ^8  [Sin  S(XG-G+Fe)  +  F#  +  XH-Hl]6  F 

+  t(-[W  +  2)(2F  +  Ge  +  8  cot  e  +  H0) -2(Ge+F) +XQF]cot  8 


+  2(F9-G)+X(X+1)G  +  XFgJsin  8 

+  ^[sin  0((Q(XF+2F  +  G0  +  G  cot  8  +  H0)  +  2(G0 +F))) ] 

+  ^[H9-H  cot  e+^-j  G0]}  9G 

{sin  8 [cot  8(He  -H  cot  8+-L  -  G0)  +  2^  F#  -H)  + 

+  x(x  +  I)H+j£?  V  +  ^  tsitl  e<H8  '«  9  +  Titi  VJ 


Integrating  by  parts  with  respect  to  0  and  0  the  respective  terms 
whose  partial  derivatives  are  broughtout,  and  using  the  negative  of 
Eq.  (A.  5): 
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JJJJsin  0{[Q(1  -X)+2](XF+2F  +  GQ  +  G  cot  6 +7J“-g  H^)  -  2x(X +2)F]6H 
+  .in  etXG-e  +  r^Hg+tjji-jF^  +  W-HjtF^ 

+  Sin  0{ [ (Q  +  2 ) (2F  +  Gg  +  G  cot  0  +  H#)  -  2(GQ +F)  +  XQF)cot  g 

-  2(Fq-G)  -X(X.  +  1)  -XFQ}ftG 

sin  0[Q(XF+2F  +  Gq+G  cot  6  +  H  )  +  2(Gq  +F) ]ft  GQ 

+  (H0-Hcote+JT^G4)6G# 

-  •‘•stcot  0(h0-h  cot  s+nihV+2(nfe  VH)+xft+l)H  + 
+  V4H  +  sin  e[H0  -B  cot  e+  TiTe  V4He 

+  tQ(XF+2F  +  G9  +  G  cot  0  +  JIJ-0  V  +  2(liTe  V°  cot  e+F)] 
wjlded, 

+  J  {{sin  0{(XG-G  +  F0)6G+[Q(XF+2F+G0+G  cot  S+^—gH^)  + 
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+  2<G6+F))6C  +  (B9  -a  cot  e  +  nH  G#)8H}ne 

+  sl"  etTlh  '7ihF«+XH-H)6r+Iih(H6-H  cot  9  + 

+  Jib  V6G  +  dbw(XF+2F+VG  «*  9*nrA> + 

+  2(~-q  H0  +  G  cot  e  +  F)]fiH}n03}ds  (A.6) 

Substituting  the  boundary  conditions,  Eqs.  (1.  lla-c),  and 
Eq.  (A-6)  into  Eq)  (1.12)  one  can  see  that  the  line  integral  of 
Eq.  (A.6)  cancels  the  line  integral  of  Eq.  (1.12).  Thus,  the  boundary 
conditions  are  automatically  included,  and  the  following  variational 
equation  results: 

JJ{{(CQ(1-X)+2](XF  +  2F+C0  +  G  cot  8+^~^H0)  -2X(\  +  2)F}6F 

+  <«-G+Fe>*Viib  (iib  VXH'H,5F(# 

+  {[(Q  +  2)(G0  +  2F  +  G  cot  e+7i~0  H  )  -2(G0+F)+XQF]cot  0 

-  2(F0-G)  -  X(X+1)G  -  \F0}6G 

+  tQ(XF+2F+Ge  +  G  cot  0+n^e%)  +2(G9  +  F)]6Ge 


+  -r  (H.  -H  cot  6  +  ■  -:1  --  G  .  )&G  . 
sin  9  8  sin  8  0'°  0 


-  t“*  e(He-Hcot  0+91^6  V+2(5ihF«'H)+x(x  +  l)H  + 


+  iih  r«)SH  +  (He  -H  cot  9+nb  C0)6He 


+  7ST  [Q(w+a  +  ee  +  B  cot  0+^  V+2(HT5  H«  + 


+  G  cot  8+F)]6H  U  sin  8  d0d0 


(A.  7) 


Note  that  Eq.  (A.  7)  does  not  have  any  second  order  derivatives 
and  that  the  notations  in  Eq.  (1.  lU)  are  those  shown  in  Eq.  (A. 7). 
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APPENDIX  B 

FINITE  ELEMENT  PROGRAM  WHEN  EIGNEVALUE,  X,  IS  REAL 
(See  Chapter  II). 
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PHOfiRAP  SINGIIAPL  >>'».  INI'UT  .OUTPUT  I  —  —  SI2.NSI  s  -  PHI  •  lull 

OlMfRSIOn  A(4>5V.A5|  ,«US9I  «A*4|65I.NGL0(I2>  .00(201  .PPOISIRI  ,AI.  I?0)  |  513. NS|  ■  -Pllll  •  !»'l 

1  tPOINI  I3»  ,Kt  I  I  Jl  .  SU.N5I  •  PHI  •  Till 

OlMtNSlON  CUHUAI6I  '  c  PAR  MAI.  DtRIVAllVE  W.H.I.  PHI  OMf.G.OR  M) 
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APPENDIX  C 

FINITE  ELEMENT  PROGRAM  WHEN  EIGENVALUE,  X,  IS  COMPLETE 
(See  Chapter  II) . 
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